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Abstract  This  work  develops  a  framework  for  SIMP-based 
topology  optimization  of  a  metallic  panel  structure  sub¬ 
jected  to  design-dependent  aerodynamic,  inertial,  elastic, 
and  thermal  loads.  Multi-physics  eigenvalue-based  design 
metrics  such  as  thermal  buckling  and  dynamic  flutter  are 
derived,  along  with  their  adjoint-based  design  derivatives. 
Locating  the  flutter  point  (Hopf-bifurcation)  in  a  precise  and 
efficient  manner  is  a  particular  challenge,  as  is  outfitting  the 
optimization  problem  with  sufficient  constraints  such  that 
the  critical  flutter  mode  does  not  switch  during  the  design 
process.  Results  are  presented  for  flutter- optimal  topologies 
of  an  unheated  panel,  thermal  buckling-optimal  topologies, 
and  flutter-optimality  of  a  heated  panel  (where  the  latter 
case  presents  a  topological  compromise  between  the  former 
two).  The  effect  of  various  constraint  boundaries,  temper¬ 
ature  gradients,  and  (for  the  flutter  of  the  heated  panel) 
thermal  load  magnitude  are  assessed.  Off-design  flutter  and 
thermal  buckling  boundaries  are  given  as  well. 

Keywords  Panel  flutter  •  Topology  optimization  • 
Aerothermoelasticity 

1  Introduction 

Panel  flutter  has  received  a  great  deal  of  attention  over 
the  last  fifty  years  within  the  computational  aeroelasticity 
community.  The  problem  is  classically  defined  by  a  very 
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thin  two-dimensional  elastic  structure,  constrained  in  some 
capacity  at  either  edge,  and  subjected  to  high  supersonic 
flow  over  the  top  surface.  For  increasing  flow  speeds,  a 
pair  of  complex  conjugate  eigenvalues  of  the  system  can 
cross  the  imaginary  axis  via  a  Hopf-bifurcation,  and  the 
subsequent  loss  of  dynamic  stability  is  defined  as  flutter. 
This  system  has  received  such  widespread  attention  due, 
in  part,  to  the  fact  that  the  aeroelastic  system  is  computa¬ 
tionally  tractable,  yet  can  display  highly-complex  coupled 
unsteady  behavior.  Furthermore,  a  wide  range  of  aerospace 
configurations  utilize  elastic  panels  (high-speed  aircraft,  jet 
engines,  space  re-entry  vehicles,  missiles,  etc.),  and  the  need 
to  prevent  aeroelastically-induced  panel  failures  can  have 
a  substantial  impact  on  the  overall  design  process.  Impor¬ 
tant  early  papers  on  this  subject  are  given  by  Fung  (1958), 
Dugundji  (1966)  and  Dowell  (1966),  more  recent  survey 
papers  are  given  by  Dowell  (1970),  Reed  et  al.  (1987),  and 
Mei  et  al.  (1999). 

Beyond  the  basic  problem  description  given  above,  many 
additional  factors  may  be  included  in  the  panel  flutter  anal¬ 
ysis  (e.g.,  a  static  pressure  differential,  various  boundary 
conditions  and  support  mechanisms,  post-flutter  nonlinear 
limit  cycle  oscillations,  acoustic  loadings:  all  of  which  are 
covered  in  the  cited  review  papers),  with  thermal  effects  of 
particular  concern  for  this  work.  If  both  ends  of  the  panel 
are  constrained  from  moving  and/or  rotating,  the  application 
of  a  temperature  field  will  weaken  the  panel,  and  a  buck¬ 
ling  temperature  can  be  identified,  above  which  the  structure 
is  statically  unstable  to  transverse  loadings.  The  heating 
fundamentally  alters  the  flutter  point  of  the  panel  as  well, 
which  is  now  defined  by  a  balance  of  inertial,  elastic,  aero¬ 
dynamic,  and  thermal  effects:  aerothermoelasticity.  Survey 
papers  on  this  general  subject  are  given  by  Garrick  (1963), 
Thornton  (1992),  McNamara  et  al.  (2008),  and  McNamara 
and  Friedmann  (2011).  Many  papers  have  been  written 
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which  specifically  detail  the  buckling  and  flutter  boundaries 
of  heated  panels:  see  Schaeffer  and  Heard  (1965),  Yang  and 
Han  (1976),  Xue  and  Mei  (1993),  and  Librescu  et  al.  (2004) 
(prescribed  temperature  field),  Abbas  et  al.  (1993)  and  Gee 
and  Sipcic  (1999)  (temperature  field  obtained  via  aerody¬ 
namic  heating),  and  Thornton  and  Dechaumphai  (1988), 
Kontinos  (1997),  Culler  and  McNamara  (2010),  and  Miller 
et  al.  (2011)  (temperature  field  obtained  via  a  balance  of 
aerodynamic  heating,  thermal  radiation,  and  transient  heat 
transfer  throughout  the  vibrating  structure). 

As  with  any  aero- structural  component,  the  weight  of 
a  panel  should  be  as  low  as  possible,  provided  that  the 
resulting  physical  response  is  acceptable.  In  the  context  of 
the  above  discussion,  this  would  imply  a  minimum-mass 
structure  with  constraints  on  buckling  and  flutter,  though 
vibrations,  deflections,  stresses,  noise,  and  fatigue  life  may 
also  be  of  concern.  Many  papers  exist  which  consider  the 
minimum  mass  of  a  thin  panel  under  a  flutter  constraint, 
or  vice-versa:  maximum  flutter  speed  under  a  mass  con¬ 
straint.  This  may  be  done  by  parameterizing  the  thickness 
distribution  throughout  the  panel  (via  either  discrete  or  con¬ 
tinuous  variables),  where  a  change  in  thickness  impacts  both 
the  inertial  and  elastic  loads  throughout  the  panel.  Early 
work  in  this  area  is  given  by  Weisshaar  (1972)  and  Peirson 
(1975),  and  additional  papers  are  found  by  Weisshaar 
(1976),  Van  Keuren  and  Eastep  (1977),  Beiner  and  Librescu 
(1983),  Barboni  et  al.  (1999),  and  Palaniappan  et  al.  (2006). 
Panel  flutter  objectives/constraints  have  also  been  handled 
with  planform  parameters  (Livne  and  Mineau  1997),  opti¬ 
mal  placement  of  piezoelectric  actuators  and  dampers  (Nam 
et  al.  1996;  Sadri  et  al.  2002;  Tanaka  et  al.  2005),  and  stack¬ 
ing  sequence  design  (Hirano  and  Todoroki  2004;  Seresta 
et  al.  2006).  Structural  optimization  of  panels  for  buck¬ 
ling  loads  is  also  well- studied  (see  Manickarajah  et  al. 
2000  and  Maalawi  2002  for  overviews  on  the  subject),  with 
thermal  effects  accounted  for  by  Foldager  et  al.  (2001), 
Chen  et  al.  (2003),  and  Wang  et  al.  (2004)  (among  many 
others).  But  the  work  of  Seresta  et  al.  (2006)  appears 
to  be  one  of  very  few  that  conducts  structural  optimiza¬ 
tion  of  a  panel  with  both  flutter  and  buckling  (thermal 
or  otherwise)  metrics,  providing  Pareto  trade-offs  between 
the  two. 

This  work  is  concerned  with  using  a  subset  of  structural 
optimization,  topology  optimization,  to  optimize  aerther- 
moelastic  panels  for  flutter  and  buckling  metrics.  Smaller 
length-to-thickness  ratios  than  might  be  traditionally  con¬ 
sidered  (L/h  ~  100)  are  utilized  here  (L/h  =  25),  so  that 
the  cross-section  of  the  two-dimensional  panel  may  be  dis¬ 
cretized  into  a  series  of  quadrilateral  finite  elements.  Each 
element  is  assigned  a  design  variable,  which  continuously 
interpolates  between  0  (void)  and  1  (solid),  and  the  opti¬ 
mal  layout  of  material  within  the  panel  may  be  obtained 


such  that  the  structure  behaves  as  intended  when  subjected 
to  thermal  and  aerodynamic  effects.  The  monograph  by 
Bendspe  and  Sigmund  (2003)  provides  a  general  overview 
of  the  method.  In  the  field  of  dynamics,  topology  opti¬ 
mization  has  been  used  to  tailor  the  response  of  structures 
subjected  to  dynamic  loads  (Ma  et  al.  1995;  Min  et  al.  1999; 
Tchemiak  2002;  Jog  2002;  Jensen  2007,  2009;  Stanford  and 
Beran  2011a,  b),  and  (of  more  concern  to  this  work)  to 
optimize  free  vibration  eigenvalues  (Pedersen  2000;  Maeda 
et  al.  2006),  frequency  gaps  (Du  and  Olhoff  2007),  or 
complete  frequency  response  spectra  (Yoon  2010a).  The 
method  has  also  been  used  for  buckling  design  (Neves  et  al. 
1995;  Rahmatalla  and  Swan  2003),  numerous  thermoelas¬ 
tic  studies  (Jog  1996;  Li  et  al.  2000;  Sigmund  2001),  and 
convection  heat  transfer  (Bruns  2007),  each  of  which  is  a 
field  that  may  be  significant  to  panel  design.  It  is  particu¬ 
larly  important  to  note  that,  given  the  large  number  of  design 
variables  utilized  in  topology  optimization  (in  comparison 
to  shape,  sizing,  or  stacking  sequence  optimization),  several 
papers  (Du  and  Olhoff  2007;  Neves  et  al.  1995)  have  noted 
that  the  large  design  freedom  can  lead  to  eigenvalue-related 
optimization  issues.  Specifically,  repeated  eigenvalues  and 
discontinuous  mode- switching  must  be  accounted  for. 

Aeroelastic  applications  of  topology  optimization  are 
rare.  Eschenauer  and  Olhoff  (2001)  and  Krog  et  al.  (2004) 
both  studied  the  topological  design  of  internal  wing  ribs, 
but  the  aerodynamic  loads  are  prescribed,  not  design- 
dependent.  Studies  which  include  aerodynamic  loading 
feedback  via  true  fluid- structure  coupling  are  given  by 
Maute  and  Allen  (2004)  (topology  of  an  underlying  wing 
structure),  Maute  and  Reich  (2006)  (compliant  airfoil  mor¬ 
phing  mechanisms),  Gomes  and  Suleman  (2008)  (wing  box 
design  for  rolling  maneuvers),  Stanford  and  Ifju  (2009) 
(membrane  wing  design),  Yoon  (2010b)  (monolithic  meth¬ 
ods  for  fluid- structure  interaction  in  topology  optimiza¬ 
tion),  Kreissl  et  al.  (2010)  (microfluiduc  devices),  and 
Stanford  and  Beran  (2011a)  (compliant  flapping  mecha¬ 
nisms).  Only  the  latter  paper  deals  with  dynamic  aeroe- 
lasticity,  rather  than  steady  fluid- structure  interaction. 
Recently,  some  researchers  have  specifically  incorporated 
aeroelastic  instability  (i.e.,  flutter)  into  the  topology  opti¬ 
mization  of  cantilevered  wing  structures:  see  Stanford  and 
Beran  (201  lb)  and  De  Leon  et  al.  (2012). 

Many  topology  and  (more  generally)  structural  opti¬ 
mization  papers  cited  above  study  the  individual  elements 
important  to  an  aerothermoelastic  panel.  To  the  best  of  the 
author’s  knowledge,  none  have  combined  them  into  a  single 
topology  optimization  framework,  which  is  the  goal  of  this 
paper.  The  development  of  such  a  framework  is  important  in 
the  sense  that  both  buckling  and  flutter  are  common  modes 
of  failure  for  aerospace  panels  (Mei  et  al.  1999),  but  design 
trends  needed  to  alleviate  them  can  be  at  strong  odds  with 
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one  another  (Seresta  et  al.  2006).  The  specific  objectives  are 

as  follows: 

1.  A  quasi-static  thermoelastic  model  will  be  developed: 
the  temperature  of  the  top  and  bottom  panel  surfaces 
is  prescribed,  and  the  conduction  through  the  topol¬ 
ogy  of  the  internal  layout  is  subsequently  computed. 
This  model  will  be  used  to  compute  stresses,  and  the 
concomitant  critical  buckling  temperatures  of  the  panel. 

2.  Assuming  that  the  characteristic  time  of  the  thermal 
system  is  large  relative  to  the  vibration  of  the  panel 
(Culler  and  McNamara  2010),  linear  piston  theory  aero¬ 
dynamics  (Ashley  and  Zartarian  1956)  (a  common  tool 
for  panel  analysis)  will  be  used  to  compute  the  flutter 
speed  of  the  panel  for  a  prescribed  top  and  bottom  panel 
surface  temperature. 

3.  The  analytical  sensitivities  of  both  the  critical  buck¬ 
ling  temperature  and  the  flutter  speed  with  respect  to 
a  large  number  of  topological  design  variables  will  be 
computed,  for  a  given  material  interpolation  scheme 
between  solid  and  void. 

4.  Three  basic  types  of  topology  optimization  problems 
will  be  solved:  maximum  thermal  buckling  load,  max¬ 
imum  unheated  flutter  speed,  and  maximum  heated 
flutter  speed  (with  a  prescribed  top  and  bottom  panel 
surface  temperature),  where  the  latter  problem  should 
pose  a  topological  compromise  between  the  two  former. 
A  set  of  constraints  will  be  carefully  posed  such  that 
mode- switching  is  avoided  during  the  optimization  pro¬ 
cess,  as  discussed  in  by  Hanoaka  and  Washizu  (1980). 

5.  The  dependence  of  the  optimal  panel  topology  upon 
various  factors  will  be  assessed:  constraint  boundaries, 
ratios  of  the  prescribed  upper  and  lower  panel  sur¬ 
face  temperatures,  and  (for  the  flutter  of  the  heated 
panel)  magnitude  of  the  thermal  load.  Furthermore, 
the  off-design  performance  of  the  topologies  will  be 
examined. 


2  Problem  description 

The  geometry  of  a  metallic  panel  structure  considered  in  this 
work  can  be  seen  in  Fig.  1.  This  test  case  is  meant  to  emu¬ 
late  the  aerospace  configurations  described  above:  an  elastic 
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Fig.  1  Panel  geometry  and  design  domain 


surface  embedded  in  an  otherwise  rigid  plane  (the  outer  sur¬ 
face  of  a  high-speed  aircraft  wing,  or  the  side  of  a  re-entry 
vehicle,  for  example),  exposed  to  flow  over  its  top  surface 
only.  The  length  and  thickness  of  the  two-dimensional  panel 
are  L  and  h ,  respectively,  the  structure  is  assumed  to  be 
infinitely-long  in  the  third  dimension,  and  both  the  lead¬ 
ing  edge  and  the  trailing  edge  are  restricted  from  moving. 
As  noted  above,  L/h  is  fixed  as  25:  this  is  smaller  than 
found  in  typical  panel  flutter  studies  (which  may  utilize 
thin  beam/plate  methods),  but  is  required  to  study  the  inter¬ 
nal  topology  of  the  structure  in  a  computationally  feasible 
manner. 

The  panel  cross  section  will  be  discretized  into  a  series 
of  bilinear  finite  elements,  and  each  element  assigned  a  den¬ 
sity  measure  x^  which  smoothly  varies  between  0  (void)  and 
1  (solid).  Portions  of  the  panel  are  fixed  as  solid  (elements 
that  lie  in  the  upper  10  %  or  the  lower  5  %,  as  drawn  in 
Fig.  1),  the  remainder  constitutes  the  design  domain.  Finite 
element  nodes  that  lie  along  the  upper  and  lower  surfaces 
are  prescribed  spatially-uniform  temperatures  T\j  and  Tl, 
respectively,  and  the  upper  surface  of  the  panel  is  subjected 
to  high-speed  supersonic  flow  (flow  density  p0 0,  flow  veloc¬ 
ity  Uoo ,  and  Mach  number  M0 0  >1).  Prescribing  the  upper 
surface  of  the  panel  to  always  remain  solid  (xe  =  1)  cir¬ 
cumvents  potential  ambiguities  concerning  the  transfer  of 
aerodynamic  loads  to  the  structure,  though  strategies  are  dis¬ 
cussed  by  Yoon  (2010b).  Prescribing  the  lower  surface  of 
the  panel  to  remain  solid  as  well,  though  this  surface  is  not 
exposed  to  aerodynamic  loads,  is  found  to  help  stabilize  the 
topology  optimization  process. 


3  Aerothermoelastic  modeling  framework 

Linear  unsteady  finite  element  solutions  are  used  to  ascer¬ 
tain  the  thermal  buckling  and  flutter  boundaries  of  the  panel 
topologies.  The  geometry  seen  in  Fig.  1  is  discretized  into 
bilinear  Q4  finite  elements,  a  common  choice  for  topology 
optimization  (Bendspe  and  Sigmund  2003).  The  remainder 
of  this  section  details  the  computational  steps  needed  to 
compute  the  relevant  metrics. 

3.1  Heat  conduction  modeling 

A  typical  assumption  in  panel  analyses  is  to  assume  a  set 
temperature  distribution  throughout  the  entirety  of  the  struc¬ 
ture  (Xue  and  Mei  1993).  This  is  reasonable  in  the  sense 
that  the  Biot  number  for  a  thin  panel  is  small  enough  to 
warrant  a  lumped-capacity  model  (Gee  and  Sipcic  1999), 
though  the  temperature  of  the  panel  could  be  linked  to 
the  flow  conditions  via  an  aerodynamic  heating  model 
(assumption  of  an  adiabatic  wall  temperature,  for  example). 
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A  compromise  between  the  two  approaches  is  taken  for  this 
work,  in  that  the  thermal  boundary  conditions  of  the  panel 
are  prescribed,  but  the  internal  conduction  is  subsequently 
computed.  Presumably,  the  small  value  of  L/h  and  the 
substantial  variations  in  mechanical  properties  throughout 
an  optimized  topology  limit  the  usefulness  of  a  lumped- 
capacity  assumption.  Furthermore,  unsteady  thermal  effects 
are  neglected  for  this  work,  assuming  that  the  time  scale  of 
the  thermal  system  is  large  relative  to  the  vibration  of  the 
panel  (Culler  and  McNamara  2010). 

In  addition  to  the  prescribed  temperature  boundary  con¬ 
ditions  stated  in  Fig.  1,  adiabatic  boundary  conditions  are 
applied  to  the  sides  of  the  panel.  The  conductivity  matrix  of 
the  panel  is  assembled  over  each  finite  element  ( e )  in  the 
typical  manner  (Cook  et  al.  2002): 


■  fe) P )  (D 

e 

where  x^  is  the  relative  density  of  each  finite  element 
and  p  is  the  penalization  power.  The  solid-void  interpo¬ 
lation  scheme  used  here  is  SIMP  (solid  isotropic  material 
with  penalization),  described  extensively  by  Bendspe  and 
Sigmund  (2003)  and  Gao  and  Zhang  (2010)  specifically 
with  regard  to  thermoelastic  problems.  Kje  is  an  elemental 
conductivity  matrix  (one  temperature  degree  of  freedom  per 
node)  for  the  fully  solid  material.  The  assembled  matrix  Kj 
is  then  partitioned  into  known  (Tc)  and  unknown  (Tx)  nodal 
temperatures: 


T  = 

'Ktxx  Ktxc~ 

pl 

u 

1  T‘  1 

where  the  prescribed  vector  Tc  is  entirely  composed  of  Tjj 
(nodes  along  the  upper  surface)  and  (nodes  along  the 
lower  surface).  For  the  linear  analysis  considered  here,  Tjj 
is  set  equal  to  unity.  The  unknowns  are  then  computed  as: 

Ktxx  •  Tx  =  Rj  (3) 

The  thermal  load  vector  Rt  is  defined  by  —Ktxc  •  Tc. 

If  Tl  is  set  equal  to  Tjj  (uniformly  heated  panel)  the  con¬ 
duction  problem  of  (3)  becomes  trivial,  as  the  temperature 
at  every  node  is  equal  to  Tjj,  entirely  independent  of  the 
topology  of  the  structure  x^.  Only  when  Tl  differs  from  Tjj 
is  thermal  conduction  activated  throughout  the  panel,  and  a 
dependency  on  x6  is  introduced. 

3.2  Panel  stresses 


initial  elastic  stresses  within  each  element,  due  to  thermal 
expansion,  are  then  given  by: 

-a  -  Te  -  E  ■  (. xe)p  - 


'  Oe 


1  -v 


1  1  0 


(4) 


where  a  is  the  coefficient  of  thermal  expansion,  v  is  the 
Poisson’s  ratio,  and  E  is  the  elastic  modulus  of  the  structure. 
A  power  law  interpolation  (SIMP)  is  used  here  as  well,  with 
the  penalty  p  assumed  to  be  identical  to  that  used  for  the 
conduction  terms  in  (1).  Furthermore,  the  thermal  expansion 
coefficient  is  assumed  to  be  independent  of  x6 :  both  of  these 
assumptions  follow  those  made  by  Sigmund  (2001).  From 
aQe ,  a  thermal  load  vector  F  may  be  computed  by  assem¬ 
bling  contributions  from  each  element.  A  linear  stiffness 
matrix  is  given  by  the  following  assembly: 


K  =  J2(Ke-  (*min  +  d  -  *min)  '  (xe)P))/(\  ~  V2)  (5) 

e 


where  Ke  is  an  elemental  stiffness  matrix  for  the  fully  solid 
material  (as  in  (1)).  This  matrix  is  computed  assuming  a 
unit  depth  of  the  structure,  and  is  divided  by  (1  —  v2)  to 
account  for  this  very-long  third  dimension  (i.e.,  a  semi¬ 
infinite  plate).  xmin  is  a  lower  bound  on  the  element  density 
(0  <  xmin  <  Xg  <  1).  In  general,  a  small  positive  value 
of  *min  is  utilized  to  prevent  singular  stiffness  matrices 
due  to  void  portions  of  the  topology,  and  also  to  allow 
for  the  reappearance,  during  the  optimization  process,  of 
material  in  elements  that  were  once  void  (i.e.,  to  prevent 
design  gradients  from  becoming  exactly  zero).  The  specific 
material  interpolation  scheme  used  in  (5)  is  a  minor  mod¬ 
ification  of  that  used  for  (1),  and  is  taken  from  Bendspe 
and  Sigmund  (2003).  The  modified  interpolation  scheme  is 
aimed  at  preventing  artificial  local  vibration  and  buckling 
modes  in  low-density  elements,  as  the  ratio  of  mass  to 
stiffness  (for  vibration  problems)  or  nonlinear  stress  stiffen¬ 
ing  to  linear  stiffness  (for  buckling  problems)  is  finite  for 
small  values  of  x^.  Additional  schemes  to  avoid  this  prob¬ 
lem  in  eigenvalue-based  topology  optimization  are  given  by 
Tchemiak  (2002),  Pedersen  (2000),  Du  and  Olhoff  (2007), 
and  Neves  et  al.  (1995).  Furthermore,  Pedersen  (2000) 
discusses  situations  where  (5)  is  unable  to  prevent  the  devel¬ 
opment  of  localized  modes,  though  the  method  is  found  to 
be  completely  satisfactory  for  the  cases  presented  here. 

Next,  a  linear  system  is  solved  for  the  thermally-induced 
deflections  of  the  structure: 


K  u  =  F  (6) 

The  stresses  produced  by  the  deflection  vector  u  are  then 
computed,  and  superposed  upon  a0e  to  obtain  the  total 
stresses  within  each  element  (Cook  et  al.  2002): 


The  temperature  is  assumed  to  be  spatially  uniform  within 
each  finite  element  (Te),  and  is  computed  by  simply  averag¬ 
ing  the  four  nodal  temperatures  found  in  the  vector  T .  The 
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where  Be  is  a  strain-displacement  matrix,  and  ue  are  the 
nodal  displacements  (two  degrees  of  freedom  per  node)  of 
each  element,  taken  as  a  subset  of  u.  It  should  be  noted  that 
a  plane-stress  formulation  is  utilized  in  (5)  and  (7),  which 
follows  from  the  common  methodology  found  in  the  panel 
literature  (Dowell  1970;  Librescu  et  al.  2004;  Culler  and 
McNamara  2010),  and  a  desire  to  reproduce  well-known 
panel  stability  boundaries  (Fig.  3).  A  plane-strain  assump¬ 
tion  is  certainly  arguable  in  light  of  the  very-long  third  panel 
dimension  (semi-infinite  plate),  though  the  boundary  condi¬ 
tions  in  this  third  dimension  are  unspecified,  both  here  and 
in  the  general  literature.  Future  work  may  assess  the  impact 
of  this  assumption  upon  the  optimal  topology. 


3.3  Thermal  buckling 


A  geometric  stress  stiffness  matrix  for  the  structure  may  be 
assembled  as: 

Ka=J2  (Ge  ■  Se  e)  '  Ge  ■  Ve)  (8) 

e 

where  Ve  is  the  volume  of  the  element,  Se  is  a  matrix  re¬ 
ordering  of  the  element  stress  ae ,  and  Ge  is  a  shape  function 
differentiation  matrix;  each  of  which  is  described  by  Cook 
et  al.  (2002).  The  buckling  eigenproblem  is  defined  as: 

(K  +  Tk.  Ka)  *k  =  0  (9) 


where  $>k  is  the  eigenvector  associated  with  the  kth  eigen¬ 
value  Tk,  and  Tk  is  a  multiplicative  factor  of  the  upper  panel 
temperature  Tjj  (Fig.  1),  which  has  been  set  to  unity  for  (2). 
Each  buckling  temperature  factor  is  positive,  and  ordered 
such  that  0  <  T\  <  T2  <  ...TNm-i  <  TNm.  Nm  is 
the  number  of  retained  modes  (for  the  purposes  of  eigen¬ 
value  separation  constraints,  as  discussed  below),  which  is 
much  smaller  than  the  total  degrees  of  freedom  of  the  elastic 
system.  The  critical  buckling  temperature  of  the  system  is 
defined  as  T*  =  T\,  which  can  be  nondimensionalized  as: 


R * 


a  ■  E  ■  T*  ■  h  ■  L2 
D 


(10) 


a  finite  strip  method,  but  their  inclusion  is  out-of- scope  for 
the  current  paper,  and  may  be  considered  in  future  efforts. 

3.4  Aerodynamic  modeling 


Linear  quasi- steady  piston  theory  relates  the  aerodynamic 
pressure  pa  at  a  point  along  the  upper  surface  of  the  panel 
to  the  instantaneous  velocity  of  that  point: 


Poo-UZo  (dw  M^-2  J_  dw\ 
y/Ml-l'  \Bx  M2,  -  l'  Uoo  '  dt  ) 


(ID 


where  x  is  the  freestream  direction  along  the  panel  in  Fig.  1, 
w  is  the  transverse  deflection  of  the  panel,  and  t  is  time. 
Terms  are  included  in  (1 1)  to  correct  the  aerodynamic  pres¬ 
sure  for  low  Mach  numbers  (Dowell  1966);  for  values  of 
Moo  much  larger  than  unity,  (11)  coincides  with  that  origi¬ 
nally  given  by  Ashley  and  Zartarian  (1956).  Equation  (11) 
has  a  stiffness-based  component  (via  dw/dx)  and  a  tempo¬ 
ral  component  (via  dw/dt ),  and  so  aerodynamic  stiffness 
and  damping  matrices  may  be  computed  by  looping  over 
each  finite  element  edge  which  lies  along  the  top  surface 
(as  only  this  surface  is  exposed  to  the  flow),  but  assembling 
into  matrices  which  are  the  same  size  as  K  and  Ka  above. 

The  panel  deflection  at  any  time  due  to  aerodynamic 
forces  is  given  by  the  vector  q ,  which  is  the  same  size  as 
the  thermal  deformation  vector  m,  and  has  two  displacement 
degrees  of  freedom  at  each  node.  The  third  and  fourth  nodes 
of  the  bilinear  finite  element  lie  along  its  upper  surface,  and 
the  third  node  lies  downstream  of  the  fourth  (i.e.,  has  a  larger 
value  of  x).  The  aerodynamic  pressure  applied  to  each  finite 
element  along  the  top  surface  can  then  be  given  by: 


Pa  = 


poo  '  Uqc 

/  —  •  [0  0  0  0  0  1  0  -1]  qe  + 


ML -2  1  1 


14-1  Uo, 


[0  0  0  0  0  1  0  l]-0, 


/ 

(12) 


where  the  bending  stiffness  of  the  panel  is  given  by  D  = 
E  -  /z3/(l2  -  (1  —  v2)).  For  a  solid  panel  with  uniform  heat¬ 
ing  (Tl  =  Tjj  =  1)  and  clamped  boundary  conditions 
(which  are  assumed  for  this  work,  as  drawn  in  Fig.  1),  an 
analytical  value  of  7?*  =  4  •  n2  is  obtained  (Shigley  and 
Mishke  2001).  Alternatively,  if  Tl  is  set  to  0  and  Tjj  =  1, 
double  this  value  is  expected:  R*  =  8  •  tt2. 

Only  buckling  modes  in  the  plane  of  the  panel  structure 
of  Fig.  1  are  included  in  the  Nm  modes,  though  buckling 
(or  flutter)  modes  can  certainly  develop  in  the  omitted  third 
direction.  These  may  be  accounted  for  within  the  context 
of  the  two-dimensional  finite  element  model  used  here  via 


where  qe  is  the  elemental  deflection  vector,  qe  is  the  time 
derivative  of  that  vector,  and  Le  is  the  length  of  the  upper 
surface  of  the  element.  It  can  be  seen  that  only  the  trans¬ 
verse  degrees  of  freedom  of  the  top  two  nodes  of  the  element 
correspond  to  an  aerodynamic  response,  in  keeping  with 
(11).  A  work-equivalent  aerodynamic  force  vector  for  that 
element  is  then  computed: 

Fae  =  -EelTi  ■  [o  o  o  o  o  1  o  l]r  (13) 

where  the  negative  sign  is  indicative  of  the  fact  that  posi¬ 
tive  aerodynamic  pressure  on  the  top  surface  will  push  the 
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structure  downward.  As  with  (11),  only  transverse  forces 
along  the  upper  surface  are  produced. 

A  global  aerodynamic  force  vector  Fa  may  be  computed 
by  assembling  over  the  finite  elements  along  the  top  surface 
(etop)«  Because  Fa  must  be  the  same  size  as  the  stiffness 
matrices  in  (9),  aerodynamic  forces  at  finite  element  nodes 
within  the  panel  are  zero.  The  force  vector  may  be  written 
as  a  linear  combination  of  q  and  q : 

Fa=J2  (F«4  =  ■  Ca  ■  q  +  X  ■  Ka  ■  q  (14) 

etop 


The  aerodynamic  stiffness  matrix  Ka  reflects  aerodynamic 
pressures  proportional  to  q  in  (12),  the  aerodynamic  damp¬ 
ing  matrix  Ca  reflects  proportionality  to  q ,  and  the  aerody¬ 
namic  pressure  parameter  is  defined  as  Dowell  (1966): 


Poo  •  Uj  -  L 3 

D • -  1 


(15) 


3.5  Panel  dynamics 

A  mass  matrix  for  the  panel  may  be  computed  as: 

M  =  J2  (Me  ■  *e)  (16) 

e 

where  Me  is  a  consistent  mass  matrix  for  the  fully  solid 
material,  and  elemental  mass  terms  are  assumed  to  be  lin¬ 
early  proportional  to  the  element  density  (Pedersen  2000). 
A  structural  damping  matrix  is  computed  via  Rayleigh 
damping  (Cook  et  al.  2002): 

C  =ac-M  +  pc-  K  (17) 

For  dynamic  compliance  (resonance)  problems,  Tcherniak 
(2002)  recommends  computing  an  elemental  Ce  based  on 
proportionality  to  Me  and  Ke,  and  then  assembling  into 
global  form  similarly  to  (16),  with  raised  to  a  power  less 
than  unity  (e.g.,  (xe)1^2).  This  can  help  remove  intermediate 
densities  from  the  topology  during  the  optimization  process, 
as  these  densities  will  dissipate  a  relatively  large  amount  of 
energy,  decreasing  actuator  quality.  This  strategy  has  been 
successfully  used  by  Stanford  and  Beran  (2011a)  as  well, 
but  for  the  eigenvalue-based  problems  considered  here  the 
simpler  approach  of  (17)  is  satisfactory. 

Combining  information  from  (5),  (8),  (14),  (16),  and 
(17),  the  equations  of  motion  for  the  panel  can  be  written  as: 

M  q  +  (C- VX  ■  Ca)  ■  q  +  (K  +  T  ■  Ka-X  ■  Ka)  •  q  =  0 

(18) 

T  is  a  multiplicative  factor  of  the  upper  panel  tempera¬ 
ture  (where  Tu ,  as  above,  is  set  to  unity).  The  applied 


temperature  T  is  entirely  unrelated  to  the  pre-computed 
critical  temperature  T*  in  (9);  values  above  or  below  T* 
may  be  selected.  Furthermore,  T  may  be  nondimension- 
alized  as  in  (10):  R  =  a  -  E  -  T  -  h  -  L2 /D.  If  both  T 
and  X  are  set  equal  to  zero,  (18)  represents  the  free  vibra¬ 
tion  of  an  unheated  panel,  and  a  characteristic  frequency 
(the  first  bending  vibration  mode)  is  identified  as  co0  = 
(4.73)2  •  4Wp  m  ’  h  -  fA),  where  pm  is  the  density  of  the 
panel.  If  a  positive  value  of  X  is  selected,  (18)  provides  the 
aerodynamically-forced  vibration  of  a  panel  with  a  mass 
ratio  of:  pi  =  p0 0  •  L/(pm  -  h).  The  inclusion  of  non-zero 
values  of  X  and  T  subjects  the  panel  to  both  aerodynamic 
and  thermal  loads,  where  it  is  understood  that  Ka  must  be 
computed  by  first  solving  for  the  temperature  distribution 
T  throughout  the  panel  (2),  then  computing  the  thermal 
deformations  u  (6),  and  finally  the  element  stresses  ae  (7). 

Defining  the  vector  rj  =  { qT  qT}T ,  (18)  may  be  placed 
in  standard  first-order  form: 


7  01  M 

[omJ  l$j 

r  o 

I 

[-K-T  ■  Ka+X-  Ka 

—C  -\-\fx  -  Ca 

(19) 


A  -  rj  =  B  -  rj 


(20) 


The  state  vector  is  assumed  to  be  of  the  form  rj  =  ifr  -  , 

which  leads  to  the  following  eigenproblem: 


(B-fo-A)-+£=  0 


(21) 


where  i/r^  is  the  right  eigenvector  associated  with  the  /cth 
eigenvalue  Pk,  both  of  which,  in  general,  will  be  complex¬ 
valued.  Specifically,  Pk  are  complex  conjugate  pairs,  where 
the  imaginary  part  is  indicative  of  the  frequency  of  vibra¬ 
tion,  and  the  real  part  dictates  stability:  a  positive  real  part 
is  an  unstable  mode  with  unbounded  growth  in  time. 


3.6  Flutter  computations 


Flutter  is  defined  as  a  loss  of  dynamic  stability  of  the  equa¬ 
tions  of  motion  (20)  about  an  equilibrium  solution,  which 
for  this  work  is  assumed  to  be  the  trivial  state:  rje  =  0. 
For  an  unheated  panel  with  a  geometry  as  drawn  in  Fig.  1, 
this  is  exactly  true,  as  the  panel  has  no  curvature  (camber) 
and  is  not  oriented  at  an  angle  of  attack  to  the  free-stream 
vector:  perturbations  to  the  system  state  will  always  set¬ 
tle  back  to  an  undeformed  state,  as  long  as  the  flow  speed 
is  below  the  flutter  speed.  For  a  heated  panel,  this  is  only 
an  assumption  (unless  both  the  heating  and  the  topology 
are  uniform),  as  there  will  be  some  static  offset  deforma¬ 
tion  associated  with  the  thermal  loads  («).  But  the  effect 
should  be  small,  and  so  this  nonlinearity  is  not  included 
(Hanoaka  and  Washizu  1980);  pre-buckling  deformations 
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have  been  similarly  ignored  in  the  thermal  buckling  eigen- 
problem  of  (9). 

Flutter  occurs  at  a  Hopf-bifurcation  point:  for  increasing 
values  of  A,  a  pair  of  complex  conjugate  eigenvalues  fa  of 
the  system  cross  the  imaginary  axis.  Defining  the  aeroelastic 
damping  of  each  mode  as  gk  =  Re  (fa),  the  system  loses 
stability  (grows  exponentially  with  time)  via  flutter  when 
gk  of  any  mode  becomes  positive  (Bisplinghoff  et  al.  1955). 
The  damping  of  the  composite  response,  G,  corresponds  to 
the  eigenvalue  with  the  largest  real  part  (the  most  unstable 
mode):  G  =  max(g&).  The  flutter  point  is  the  point  where 
G  =  0,  and  occurs  at  A,  =  A,*;  the  imaginary  portion  of 
this  critical  eigenvalue  at  A,*  is  the  flutter  frequency,  <u*.  An 
example  of  this  process  is  given  in  Fig.  2.  The  eigenvalues  at 
A.  =  0  (wind  off)  are  the  uncoupled  vibration  eigenvalues  of 
the  system,  and  Re  (fa)  is  entirely  governed  by  the  structural 
damping  C.  If  both  damping  parameters  (ac  and  fa)  are  set 
to  zero  (which  is  not  the  case  in  Fig.  2),  G  =  0  at  this  point 
as  well,  though  A.  =  0  is  not  a  Hopf-bifurcation  point.  If  the 
applied  temperature  T  is  larger  than  the  critical  value  of  T*, 
the  buckled  panel  will  have  at  least  one  eigenvalue  whose 
imaginary  portion  is  equal  to  zero  at  A.  =  0. 

Several  methods  can  be  used  for  locating  the  flutter 
point  in  a  direct  manner:  most  require  an  initial  guess  for 
A.*,  <z>*,  and  the  eigenvector  \/rR*  (aeroelastic  mode)  asso¬ 
ciated  with  this  critical  eigenvalue.  This  can  be  done  by 
setting  A.  =  0,  computing  the  eigenvalues  fa ,  and  evaluating 
G  =  ma x(Re(fa)).  If  G  is  less  than  0,  A.  is  augmented  by 
some  A  A.,  and  the  process  is  repeated  until  G  becomes  pos¬ 
itive.  Provided  that  A  A.  is  not  too  large,  the  value  of  A.  at 
which  this  occurs  should  be  a  reasonable  approximation  to 
A.*,  and  the  least-stable  mode  provides  initial  guesses  for 
co *  and  i/rR*.  The  repetitive  eigenvalue  computations  dur¬ 
ing  the  search  for  A.*  can  be  very  costly,  particularly  if  A  A. 
is  of  a  small  size.  In  an  optimization  context,  the  process 
may  be  conducted  just  once  for  the  initial  baseline  design: 
a  subsequent  design  iterate  may  use  the  previous  design 
iterate’ s  flutter  point  as  an  initial  guess.  It  will  be  shown 
however,  that  eigenvalue  information  is  needed  in  the  range 
0  <  A.  <  A.*  to  properly  stabilize  the  flutter  optimization 
problem,  and  so  this  A, -marching  process  is  conducted  for 
every  design  iterate. 


Once  an  initial  guess  has  been  obtained,  the  precise  loca¬ 
tion  of  the  flutter  point  may  be  computed  with  a  class  of 
techniques  that  track  the  position  of  the  least  stable  eigen¬ 
mode.  Three  methods  fall  within  this  class:  (1)  direct  eige- 
nanalysis  for  the  most  unstable  mode,  (2)  direct  analysis  of 
the  expanded  system  of  equations  for  the  Hopf-bifurcation 
point,  and  (3)  an  inverse  power  method  with  shifting. 
Only  the  first  will  be  discussed  here  due  to  its  inher¬ 
ent  simplicity;  the  interested  reader  may  study  the  second 
approach  (Griewank  and  Reddien  1983;  Morton  and  Beran 
1999),  or  the  third  approach  (Timme  et  al.  2011),  which 
have  both  been  successfully  applied  to  complex  aeroelastic 
problems. 

The  first  approach  uses  Newton’s  recurrence  formula  to 
drive  G  to  zero: 

Xn+x  =  Xn  —  To  ■  (9G/3X)-1  •  G  (A")  (22) 

where  n  is  the  iteration  number.  Solving  the  eigenproblem 
(21)  at  the  current  iterate  Xn  provides  G(A7),  which  is  the 
damping  of  the  least  stable  mode  (if  the  A, -marching  process 
detailed  above  is  used,  G(A.°)  will  be  positive).  The  right 
and  left  eigenvectors  may  also  be  obtained  (irR,  i?L),  and 
then  the  slope  of  the  damping  with  respect  to  A.  is  computed 
as  (Murthy  and  Haftka  1988): 


where  it  can  be  seen  in  (20)  that  A  has  no  dependence  on 
A.  (dA/dX  =  0),  and  the  mode  number  k  in  (23)  corre¬ 
sponds  to  the  least  stable  mode.  The  region  of  convergence 
of  (22)  is  wide  above  the  flutter  point  (Xn  >  A,*,  where  the 
unstable  mode  is  clearly  distinguished  from  the  others,  as 
seen  in  Fig.  2),  but  narrow  below  A,*.  As  such,  an  under¬ 
relaxation  factor  oo  is  typically  needed  to  prevent  overshoots 
that  may  drive  the  approximation  well  below  the  flutter 


Fig.  2  Normalized  eigenvalue 
migration  for  a  solid  panel  with 
a  volume  fraction  of  0.35 
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point.  This  process  requires  the  repeated  computation  of  a 
right  and  a  left  eigenvector,  whose  cost  will  depend  strongly 
upon  the  sparsity  of  A  and  B  (Lehoucq  et  al.  1998).  Upon 
convergence  of  (22):  X  X*,  go  of,  i/rR  i/rR*, 
and  \/rL  ^rL*.  Information  about  higher  eigenmodes  at 
the  flutter  point  (/3^,  i/rR*,  ^*)  are  not  directly  available 
using  this  method,  but  they  can  be  obtained  by  re-solving 
(21)  at  A*. 

A  verification  study  is  given  in  Fig.  3,  where  data  from 
the  aeroelastic  model  formulated  above  is  compared  to  data 
from  Erickson  (1966)  and  Beloiu  et  al.  (2005).  This  is 
done  for  a  solid  panel:  referencing  Fig.  1,  each  element  has 
a  density  of  unity.  Furthermore,  the  ratio  of  the  mass 
ratio  to  the  Mach  number,  /x/Moo,  is  set  to  0.1  (a  typical 
value  for  aeronautical  applications  (Dowell  1966)),  and  the 
structural  damping  C  is  set  to  zero.  Two  stability  bound¬ 
aries  are  shown  in  Fig.  3.  The  upper  boundary  shows  the 
(nearly-linear)  decrease  in  the  flutter  speed  (X*)  as  heat¬ 
ing  is  uniformly  applied  to  the  structure  (i.e.,  increased  R , 
the  normalized  value  of  the  multiplicative  factor  T  in  (18)). 
The  edges  of  the  panel  are  constrained  from  moving,  and  so 
the  compressive  thermal  stresses  weaken  the  structure  and 
thus  lower  its  flutter  speed.  The  lower  set  of  points  in  the 
figure  details  the  boundary  across  which  a  buckled  panel, 
for  a  given  value  of  R ,  is  blown  flat  by  the  aerodynamic 
loads.  This  boundary  emanates  from  R  =  4  •  tt2,  noted 
above  as  being  the  critical  buckling  load  7?*  of  the  panel 
(Shigley  and  Mishke  2001).  Except  for  7?*,  this  latter 
boundary  is  of  no  interest  for  the  current  topological  design 
study,  but  does  present  a  useful  comparison  with  published 
data. 

All  combinations  of  (. X ,  R)  below  the  flutter  boundary 
and  above  the  buckling  boundary  would  ultimately  result  in 
a  flat  panel.  Data  points  below  the  lower  boundary  result 
in  a  static  buckled  panel,  and  data  points  above  the  upper 
boundary  result  in  nonlinear  limit  cycle  oscillations  (which 
may  or  may  not  be  periodic,  depending  on  whether  R  is  to 


Fig.  3  Flutter  and  buckling  boundaries  for  a  solid  panel  with  a  vol¬ 
ume  fraction  of  1,  as  compared  with  Erickson  (1966)  and  Beloiu  et  al. 
(2005) 


the  left  or  right  of  the  point  at  which  the  two  boundaries  in 
Fig.  3  meet)  0,  but  this  behavior  is  beyond  the  scope  of  the 
current  work.  Erickson  (1966)  and  Beloiu  et  al.  (2005)  use 
Euler-Bernoulli  thin-beam  modeling,  with  clamped  bound¬ 
ary  conditions  at  either  end,  to  compute  these  boundaries, 
whereas  the  current  work  uses  a  lattice  of  plane-stress  bilin¬ 
ear  finite  elements.  Despite  this  difference,  and  despite  the 
fact  that  L/h  seen  in  Fig.  1  is  clearly  too  large  to  be  consid¬ 
ered  a  thin  beam/panel  (much  larger  values  are  presumably 
used  in  the  two  references),  the  current  model  matches  well 
with  the  published  results.  Unheated  flutter  speeds  are  given 
as  636.57  by  Erickson  (1966),  and  computed  as  633.05  in 
the  current  work.  The  buckling  load  R* /tt2  is  less  accurate, 
computed  as  3.906  for  the  current  work,  whereas  Beloiu 
et  al.  (2005)  arrive  closer  to  the  analytical  solution  of  4. 
Numerical  experiments  have  indicated  that  this  is  not  due  to 
discretization  errors,  but  rather  to  the  geometrical  violation 
of  the  thin  beam  assumption.  This  same  problem  leads  to 
a  moderate  under-prediction  of  the  flutter  speed  for  heated 
panels. 

4  Sensitivity  analysis 

Topological  design  studies  formulated  below  utilize  the 
thermal  buckling  eigenvalues  7&,  the  flutter  speed  X*,  and 
the  dynamic  eigenvalues  fik  for  flow  speeds  in  the  range 
0  <  X  <  X*,  to  form  various  objective  functions  and 
constraints.  The  number  of  design  variables  (x,  a  vector 
which  collects  the  design  density  for  each  finite  ele¬ 
ment)  will  be  very  large,  and  so  design  derivatives  must  be 
computed  analytically,  as  is  typically  the  case  for  topology 
optimization. 

4.1  Buckling  derivatives 

Considering  the  buckling  eigenvalue  problem  of  (9),  the  lin¬ 
ear  stiffness  matrix  K  is  an  explicit  function  of  the  element 
densities  x ,  though  the  (nonlinear)  stress-stiffness  matrix 
Ka  is  an  implicit  function  of  the  nodal  temperatures  and 
element  stresses  as  well,  and  of  course  these  dependen¬ 
cies  must  be  accounted  for.  The  total  derivative  of  this 
matrix  is: 

dKa  =  dKa  dKa  du  dKa  dTx 

dx  dx  du  dx  d  Tx  dx 

Explicit  x  -dependencies  come  from  both  a  0e  and  the  second 
term  in  (7).  Implicit  dependencies  via  the  thermal  deforma¬ 
tion  u  are  due  to  the  second  term  in  (7),  and  dependencies 
via  the  unknown  nodal  temperatures  Tx  are  due  to  the 
aQe  computation  in  (4)  (where  a  linear  interpolation  from 
nodal  values  Tx  to  elemental  values  Te  is  assumed  in  the 
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chain  rule  computation  of  3 Ka /3 Tx).  The  derivatives  of  (3) 
and  (6)  are: 


9KTxi 

dx 


Tx  +  Kj 


dTx 

dx 


3  Rj 
dx 


(25) 


of  repeated  eigenvalues  may  be  computed  (Seyranian  et  al. 
1994),  but  the  current  work  instead  forces  adjacent  eigen¬ 
values  to  be  separated  by  a  certain  amount  during  the 
optimization  process.  Additional  discussion  regarding  this 
point  is  provided  below. 


3  K  3  u  dF  3  F  3  F  dTx 

- u  +  K  —  =  —  =  —  + - -  (26) 

dx  dx  dx  dx  3  Tx  dx 

Explicit  derivatives  of  the  thermal  load  vector  Rj  are 
computed  via  the  matrix  K  jxc .  Derivatives  of  F  may  be 
computed  via  the  a  Qe  computation  of  (4). 

With  this  information,  the  gradients  of  the  buckling 
eigenvalues  are: 


dKa  m  du_  .  dKa  <  dT 2 
d U  dx  '  d Tx  dx 


dTk  _  (®k)T  ■  (ff"  +  Tk  ■  (jf  + 

(**)r- «*•(♦*) 

T  ( 8K  du  dF  dF  dT 

+  (K)  - u  +  K - - 

V  dx  dx  dx  dT 


))•(**) 


+  (»r)' 


—^■Tx  +  Kz  ■  — 

dx 


dx 
3  Rt 


dx  dx 


(27) 


The  first  term  in  (27)  is  the  well-known  derivative  of 
an  eigenvalue  (Murthy  and  Haftka  1988)  (where  (4>&)r  • 
Ka  -  (Q>k)  is  a  modal  normalization  condition),  the  second 
two  represent  the  multidisciplinary  nature  of  the  problem 
(Sigmund  2001):  Xu  and  Xj  are  adjoint  vectors,  both  of 
which  pre-multiply  terms  which  are  exactly  equal  to  zero. 
Collecting  terms  that  pre-multiply  du/dx  provides  a  linear 
system  of  equations  for  Xu : 


K-Xu  =  -Tk 


c**)r  •  •  (**) 

(<S>k)T  ■  Ka  ■  (**) 


(28) 


Similarly,  Xj  is  computed  by  collecting  terms  associated 
with  dTx/dx\ 


4.2  Flutter  derivatives 

Gradients  of  the  dynamic  eigenvalue  problem  (21)  pro¬ 
ceeds  in  a  similar  manner  to  the  process  given  above, 
though  differences  arise  because  the  dynamic  eigenprob- 
lem  is  complex- valued  (i.e.,  left  and  right  eigenvectors  are 
needed),  and  also  due  to  dependence  on  the  dynamic  pres¬ 
sure  parameter  X.  If  eigenvalue  derivatives  are  desired  at  a 
given  value  of  X  (less  than  the  flutter  speed  X*),  these  may 
be  computed  as  Murthy  and  Haftka  (1988): 


dPk 

dx 


■  ( 


dB  ■  dB  du  ■  dB  dTx 
3jc  '  du  ’  dx  '  d  Tx  '  dx 


W  ■*■(*?) 


#  /3 K  du  dF  dF  dT 

T  (Vu)  •  (  - ■  u  ~b  K  - - -  — 

V  J  1  3jc  dx  dx  dT 


dx 


+  (Vt) 


dx 


•Tx  +  KTja 


dT  x  dR7 


dx  dx 


(31) 


where  vu  and  vT  are  adjoint  vectors.  It  can  be  seen  that  A, 
which  is  only  composed  of  the  mass  matrix,  only  has  an 
explicit  dependence  on  x ,  whereas  B  depends  on  the  ther¬ 
moelastic  problem  as  well,  due  to  the  inclusion  of  K a  in  this 
matrix  ((19)  and  (20)).  Similar  to  (28)  and  (29),  the  adjoint 
vectors  are  given  by: 


K  vu  =  — 


i_jj :  m 

T  -  A  •  (t?) 


(32) 


K  7 


dF 

dfl 


(*k)T  ■  •  (**) 

■K- Tk •  ^  dL*  -  (29) 

(®k)T  ■  Ka  ■  ($*) 


Having  solved  these  two  systems  (once  per  each  eigenvalue 
Tk  of  interest),  (27)  can  be  written  in  final  form: 


dTk 

dx 


T  +  Tk-*fc)-(*k) 


+  (ItY 


T  -Ka- 

■(*0 

dK 

dF 

dx 

dx 

9Ktxx 

Tx- 

dx 


dRj 

dx 


(30) 


KTxx  •  vT  = 


dF 

~dT~x 


Vu  ~ 


dT  x 


o rt )  A 

and  the  final  eigenvalue  derivatives  are: 

m  (+f  -  (S  -  ft  -  If)  ■  (*?) 

dx 


r)  -a ■  (t?) 


,  ,'dK  dF 

+  (vu)  ■  -T" 

dx  dx 


+  (VtY 


d_K^ 

dx 


•  Tx  — 


dRj 

dx 


(33) 


(34) 


It  should  be  noted  that  (30)  specifically  assumes  that  none 
of  the  eigenvalues  Tk  are  repeated  (i.e.,  bimodal).  Gradients 


As  above,  this  formulation  assumes  that  none  of  the  eigen¬ 
values  Pk  are  repeated. 
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If  eigenvalue  derivatives  are  desired  at  the  flutter  speed 
A*,  then  a  A-dependence  must  be  included  in  the  chain  rule: 

dJl  =  m  dJl  ar 

dx  dx  3  A.  dx 

where  the  partial  derivative  df3£/dx  may  be  computed  using 
(34),  with  all  terms  evaluated  at  A*.  The  derivative  of  the 
eigenvalues  with  respect  to  A  has  already  been  given  in  (23), 
but  will  be  provided  in  a  more  general  form  here: 

M  =  (^L*)r  •  Mr  •  (ID  (36) 

dx  (+k*)T  ■*■(*?*) 

The  flutter  point  derivative  dX*/dx  is  required  because 
of  (35),  but  also  because  A*  will  be  used  as  a  topologi¬ 
cal  objective  function.  It  may  be  computed  by  considering 
the  aeroelastic  damping  at  the  flutter  point,  which  is  by 
definition  zero,  regardless  of  the  design  vector: 


G*  =  max  (Re  (/$£))  =  0 


dG* 

dG*  dG 

dx*  „ 

dx 

- =  0 

(37) 

dx  dX 

dx 

d  A* 

(  3G\_1 

dG* 

=  -  — 

(38) 

dx 

\dx 

dx 

The  term  3G/3A  is  the  real  part  of  (36)  (which  has  already 
been  computed  for  use  in  (22)  during  the  flutter  point 
search),  and  the  term  dG* /dx  is  the  real  part  of  (35),  with 
the  mode  number  k  in  both  equations  corresponding  to  the 
flutter  mode. 


5  Optimization  strategies 


Of  the  various  restriction  methods  used  in  topology  opti¬ 
mization,  a  spatial  filter  of  the  design  gradients  is  the  most 
popular  (Bendspe  and  Sigmund  2003).  For  the  heavily- 
constrained  optimization  problems  developed  below,  this 
degradation  in  gradient  accuracy  was  found  to  cause  dif¬ 
ficulties  in  improving  the  objective  function  in  a  feasible 
manner.  As  such,  a  density  filter  developed  by  Bruns  and 
Tortorelli  (2001)  is  used  for  this  work,  where  the  element 
densities  x  are  obtained  by  filtering  a  vector  of  design 
variables  x.  This  is  done  via  a  linearly-decaying  cone- 
shape  function,  generically  represented  by  a  filtering  matrix: 
x  =  H,  -  x.  The  gradients  may  then  be  altered  in  an  exact 
fashion  via  the  chain  rule.  For  example: 


3A*  _  3A* 
dx  dx 


(39) 


where  the  gradients  with  respect  to  x  are  what  is  com¬ 
puted  in  the  previous  section,  but  gradients  with  respect  to 


the  actual  design  variables  x  are  required  by  the  optimizer. 
A  continuation  method  is  then  used:  once  the  optimization 
problem  has  converged  to  within  a  suitable  tolerance,  the  fil¬ 
tering  radius  is  decreased,  and  the  problem  is  restarted.  As 
the  radius  becomes  zero,  x  ^  x.  This  was  found  to  be  a 
suitable  strategy  for  obtaining  a  final  topology  largely  com¬ 
posed  of  strictly  solid  and  void  elements.  The  initial  radius 
size  is  set  to  15  %  of  the  panel  thickness  h. 

The  first  topology  optimization  problem  considered  here 
seeks  to  maximize  the  critical  buckling  temperature: 

max 
x 

0  <  *min  <  Xe  <  l 

s.t.  :  vT  x<V* 

Tk  ~  Tk-i  >  st 

where  the  elemental  design  variables  are  constrained  to  lie 
between  x min  and  unity.  Because  the  filter  H  is  volume¬ 
preserving,  the  elemental  densities  x  will  also  lie  within 
these  side  constraints.  The  volume  constraint  V*  is  required 
as  an  implicit  penalty  on  intermediate  densities  via  the 
SIMP  scheme  (Bendspe  and  Sigmund  2003),  and  is  writ¬ 
ten  directly  in  terms  of  x  rather  than  x.  The  final  set  of 
constraints  in  (40)  requires  that  the  first  Nm  eigenvalues  be 
separated  by  some  small  tolerance  sj,  in  order  to  prevent 
discontinuities  associated  with  critical  mode  switching.  For 
the  high  aspect  ratio  of  the  panel  cross  section  in  Fig.  1,  this 
constraint  is  typically  active  only  for  higher  mode  numbers 
k  (as  opposed  to  the  closely  spaced  first  and  second  buckling 
eigenvalues  seen  during  optimization  of  lower  aspect  ratios 
by  Neves  et  al.  1995),  and  is  of  only  minor  importance. 

The  second  topology  optimization  seeks  to  maximize  the 
flutter  speed  A*.  As  seen  in  Fig.  2,  flutter  may  occur  (as  the 
parameter  A  is  increased)  when  the  imaginary  portions  of 
two  distinct  eigenvalues  coalesce.  Shortly  after  this  coales¬ 
cence,  the  real  part  of  one  of  these  eigenvalues  (mode  1, 
in  the  figure)  becomes  positive,  defining  A*.  Simply  max¬ 
imizing  A*  during  the  optimization  process  will  inevitably 
cause  critical  mode  switching,  where  the  eigenvalues  of  two 
higher  modes  coalesce,  redefining  the  flutter  point.  This 
may  happen  at  a  value  of  A  well  below  the  previous  design 
iterate’s  flutter  point  A*,  causing  a  sharp  drop  in  the  objec¬ 
tive  function  (Odaka  and  Furuya  2005).  In  order  to  prevent 
this,  a  series  of  constraints  can  be  formulated  such  that  the 
imaginary  portions  of  two  consecutive  eigenvalues,  fik-l 
and  Pk,  be  separated  by  some  finite  amount  at  every  value 
of  A  between  0  (no  flow)  and  A*  (Langthjem  and  Sugiyama 
1999): 


e  =  l 


(40) 


k  =  2,...,Nm 


A  cok=  min  Im  (fa  (A)  —  fa- 1  (A))  >  s^  k  =  3,...,Nm 

0<k<k* 

(41) 
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where  the  dynamic  eigenvalues  have  been  ordered  based 
upon  their  imaginary  parts.  The  point  of  minimum  separa¬ 
tion  between  two  modes  is  XAoJk,  and  may  coincide  with  A.*. 
If  the  A. -marching  method  described  above  is  used  to  locate 
an  initial  guess  for  the  flutter  point,  XA(°k  and  A  cok  for  each 
mode  pair  will  be  available  at  the  end  of  this  process.  These 
will  only  be  approximate  values,  but  provided  A  A.  is  of  a 
moderate  size,  their  accuracy  will  be  acceptable  for  the  cur¬ 
rent  purpose,  and  no  direct  method  to  hone  in  on  an  exact 
value  of  XA(°k  is  considered  here. 

The  frequency  separation  A&>2  is  of  no  consequence  (as 
interactions  between  modes  1  and  2  will  always  cause  flut¬ 
ter,  and  XAc°2  is  always  equal  to  A.*),  and  so  k  starts  at  3  in 

(41) .  If,  for  two  consecutive  eigenvalues,  this  frequency  sep¬ 
aration  is  minimum  at  the  flutter  point  (XAa)k  =  A.*),  then 
the  constraint  gradients  may  be  computed  with  (35).  If,  on 
the  other  hand,  the  point  of  minimum  separation  occurs  at 
any  other  point  (XAc°k  <  A,*),  the  second  term  in  the  chain 
rule  of  (35)  becomes  zero.  This  is  because  (41)  is  a  local 
minimum:  d(A(Ok)/dX  =  0.  Similar  arguments  can  be  made 
(Greene  and  Haftka  1991)  for  critical  point  constraints  in 
time  during  a  transient  analysis. 

For  variable-thickness  aeroelastic  plate  problems,  Odaka 
and  Furuya  (2005)  show  that  the  frequency- separation  con¬ 
straint  becomes  more  important  for  larger  side-bounds  on 
the  thickness  of  each  finite  element.  Topological  design 
represents  an  extreme  of  this  process  (as  elements  can 
become  completely  void),  and  so  this  set  of  constraints  is 
expected  to  be  very  important.  If  flutter  is  always  char¬ 
acterized  by  a  modal  coalescence,  then  the  constraints 
of  (41)  are  enough  to  stabilize  the  optimization  problem, 
as  the  “nature”  of  the  eigenvalue  migration  is  forced  to 
remain  unchanged  during  the  design  process  (Langthjem 
and  Sugiyama  1999).  For  the  topologically-optimized  pan¬ 
els  considered  here,  however,  it  was  found  that  flutter  could 
occasionally  occur  without  a  strong  eigenvalue  coalescence. 
For  these  cases,  the  enforcement  of  a  series  of  A cok  con¬ 
straints  is  not  enough  to  prevent  mode  switching  during  the 
optimization  process,  as  the  interaction  between  modes  k 
and  k  —  1  may  still  cause  flutter  even  though  A  cok  >  £<y. 
Additional  constraints  are  then  required  to  stabilize  the 
process: 

Re  (ft*)  =  g*k  <  gk\x=0  k  =  2,...,Nm  (42) 

where  it  is  required  that  the  aeroelastic  damping  at  the  flutter 
point  be  more  stabilizing  (i.e.,  more  negative)  than  the  struc¬ 
tural  (no  flow)  damping  (similar  constraints  are  considered 
by  Haftka  1975).  This  constraint  is  impossible  to  satisfy  for 
mode  1,  as  g*  =  0  is  the  definition  of  a  flutter  point,  and  so 

(42)  starts  with  mode  2. 

The  combination  of  frequency  separation  (41)  and  damp¬ 
ing  separation  (42)  was  found  to  prevent  mode  switching 


for  all  of  the  cases  presented  in  this  work.  The  complete 
optimization  problem  is  written  as: 


max^* 

x 
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II 
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Both  the  thermal  buckling  and  the  flutter  optimization 
problems  are  solved  with  the  method  of  moving  asymp¬ 
totes  (Svanberg  1987),  a  scheme  known  to  be  effective 
for  topology  problems  with  large  numbers  of  variables  and 
constraints. 

Presumably,  it  may  be  expected  that  an  alternative  han¬ 
dling  of  the  eigenvalue- switching  issues,  as  opposed  to 
the  separation  constraints  in  (40)  and  (43),  could  pro¬ 
vide  superior  optimal  results  (larger  values  of  T*  or  A.*) 
by  removing  limitations  on  the  feasible  design  space.  For 
the  thermal  buckling  problem  of  (40),  a  bound  formula¬ 
tion,  in  conjunction  with  an  eigen-tracking  scheme  (when 
the  modes  inevitably  switch)  and  algorithms  to  compute 
bimodal  eigenvalue  sensitivities  (Seyranian  et  al.  1994)  may 
be  a  preferred  method.  As  noted  above,  however  (and  shown 
in  Fig.  12),  switching  of  the  critical  thermal  buckling  mode 
is  never  observed  for  this  particular  panel  geometry,  and 
the  separation-based  formulation  in  (40)  should  give  very 
similar  results  to  a  bound  formulation. 

The  separation  constraints  of  the  aeroelastic  optimiza¬ 
tion  problem  in  (43),  on  the  other  hand,  will  be  shown  to 
be  highly  active,  and  unquestionably  drive  the  design  pro¬ 
cess  to  some  extent.  However,  alternative  schemes  to  those 
used  in  (41)  and  (42)  are  less  obvious,  due  primarily  to  the 
fact  that  the  flutter  problem  is  an  eigenvalue  migration  prob¬ 
lem,  rather  than  just  computing  and  optimizing  eigenvalues 
at  some  fixed  point.  Mode  tracking  schemes  and  repeated 
eigenvalue  sensitivity  methods  remain  viable,  but  a  bound 
formulation  becomes  far  more  complex.  Multiple  flutter 
points  would  be  obtained  to  formulate  the  series  of  bound 
constraints:  A.*  <  A,|  <  .  ..A/j^  j  <  A.^.  The  computa¬ 
tional  cost  of  this  is  large,  as  the  range  of  the  A. -marching 
method  (to  obtain  initial  guesses  for  each  A.p  would  need 
to  extend  to  very  high  flight  speeds,  and  then  the  recurrence 
formula  of  (22)  used  to  strongly  converge  to  each  flutter 
point.  It  is  also  very  possible,  during  the  optimization  pro¬ 
cess,  for  the  slope  dG/dX  of  a  higher  mode  to  become  zero 
at  the  flutter  point,  and  that  mode  then  cease  to  flutter  at  all 
(a  hump  mode,  Haftka  1975).  This  would  pose  difficulties 
for  the  bound  method,  in  the  sense  that  the  number  of  flutter 
points  Nm  would  change. 
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The  separation  constraints  of  the  current  flutter  opti¬ 
mization  method  are  likely  to  limit  the  design  process  to 
some  extent,  but  only  require  the  computation  of  a  single 
(lowest)  flutter  point,  and  are  therefore  relatively  inex¬ 
pensive.  Clearly,  more  research  is  needed  in  this  area  to 
compare  these  methods,  not  just  for  topology  optimization, 
but  for  any  design  of  fluttering  structures  that  relies  upon 
gradients-based  optimization. 


6  Results 

All  results  presented  here  are  for  a  titanium  panel  (a  typical 
material  choice  for  thermoelastic  panel  flutter  studies,  see 
Librescu  et  al.  2004  and  Culler  and  McNamara  2010)  oper¬ 
ating  at  fi/M0 o  =0.1  and  an  aspect  ratio  L/h  =  25.  The 
damping  parameters  are  set  as  ac  =  5  s-1  and  =  10-7  s. 
The  penalization  factor  p  is  set  to  3,  the  first  8  eigenmodes 
( Nm )  are  considered  for  the  separation  constraints  in  (40) 
and  (43),  and  xmin  is  set  to  0.01.  Three  basic  sets  of  results 
are  given:  flutter  of  an  unheated  panel,  buckling  of  a  heated 
panel,  and  flutter  of  a  heated  panel,  where  the  latter  prob¬ 
lem  should  pose  a  topological  compromise  between  the  two 
former.  For  each  optimal  topology,  a  logical  comparison  is 
to  an  entirely  solid  panel  with  the  same  mass.  This  is  drawn 
in  Fig.  4  for  a  given  volume  fraction  V*,  where  finite  ele¬ 
ments  in  the  upper  portion  of  the  panel  are  set  as  solid;  the 
remainder  are  void.  The  domain  thickness  h  is  still  used  as  a 
length-scale  for  normalization  in  (10)  and  (15)  (rather  than 
V*  •  h),  in  order  to  maintain  a  consistent  comparison  with 
the  optimal  topologies. 

6.1  Flutter  of  an  unheated  panel 

Results  are  first  obtained  for  a  volume  fraction  V*  of  0.35, 
and  a  frequency  separation  constraint  (e^,  normalized  by  the 
characteristic  frequency  co0)  of  0.95.  The  panel  is  unheated, 
a  situation  which  is  obtained  by  setting  the  multiplicative 
temperature  T  to  zero  in  (18).  The  topology  of  the  panel  at 
selected  design  iterations  of  (43)  is  given  in  Fig.  5,  while 
the  actual  objective  function  and  constraint  metrics  are  plot¬ 
ted  at  every  design  iteration  in  Fig.  6.  The  baseline  design 
in  Fig.  5  consists  of  setting  the  density  of  each  element  in 
the  design  domain  to  a  value  slightly  less  than  0.35:  account¬ 
ing  for  the  areas  of  the  panel  which  are  fixed  as  solid  (xe  =  1 


solid 


l 


Moo;  U*,;  pm 


V*-h 


void 


Fig.  4  Geometry  of  a  fully-solid  panel  with  a  volume  fraction  of  V* 


within  the  top  10  %  and  bottom  5  %  of  the  panel),  the  vol¬ 
ume  fraction  of  the  complete  structure  vT  •  x  is  then  exactly 
equal  to  0.35.  Within  the  first  30  iterations,  material  is  allo¬ 
cated  to  the  center  and  the  two  ends  of  the  panel,  where  the 
former  is  mostly  due  to  inertial  considerations,  and  the  latter 
due  to  stiffness  effects.  The  structure  is  entirely  symmetric 
at  this  point,  though  the  biasing  effects  of  the  flow  vector 
(plotted  at  the  top  of  Fig.  5)  provides  the  optimizer  with 
some  incentive  to  develop  asymmetric  topologies  in  order 
to  further  increase  the  flutter  point  X*.  The  advent  of  asym¬ 
metry  after  30  iterations  occurs  (not  coincidentally)  at  the 
same  point  when  two  of  the  frequency  separation  constraints 
(A&>5  and  A<z>6)  become  active,  and  convergence  becomes 
much  slower  as  the  optimizer  must  travel  along  these  highly 
nonlinear  constraint  boundaries. 

A  third  frequency  separation  constraint  (Aco^)  becomes 
active  at  iteration  380,  by  which  time  the  skewed  shape  of 
the  center  mass  and  the  cross-bracing  at  the  leading  edge 
is  fully  formed.  The  topology  at  the  trailing  edge  contin¬ 
ues  to  develop  through  to  iteration  700,  at  which  point  the 
entire  process  has  converged.  Constraints  A&q,  Aco^,  and 
A<z>8  (the  latter  is  very  large,  and  not  shown  in  Fig.  6) 
remain  inactive,  as  do  all  of  the  damping  separation  con¬ 
straints  (42).  For  this  case,  the  frequency  separation  has 
acted  as  a  more-conservative  surrogate  for  the  damping  sep¬ 
aration  in  preventing  critical  mode  switching,  but  for  many 
of  the  cases  presented  below,  both  sets  of  constraints  are 
active.  Though  not  shown  in  Fig.  6,  the  volume  constraint 
is  also  active  beyond  iteration  20.  After  iteration  700,  the 
continuation  process  described  above  is  utilized,  by  contin¬ 
ually  shrinking  the  radius  of  the  filter  H  and  re-running  the 
optimizer  to  convergence.  Only  the  final  design  from  this 
process  (when  the  filter  radius  is  zero)  is  given  at  the  bottom 
of  Fig.  5:  the  structure  has  largely  converged  to  a  solid- void 
design  by  removing  the  blurred  boundaries  of  the  topology, 
as  intended.  The  flutter  speed  has  been  further  increased 
during  this  process,  and  the  active  constraints  at  iteration 
700  remain  active.  All  topologies  shown  in  Fig.  5  are  the  fil¬ 
tered  element  densities  jc,  as  opposed  to  the  design  variables 
x ,  which  have  no  physical  meaning.  For  the  final  design 
obtained  via  continuation,  the  filter  radius  is  zero,  and  so 
x  and  x  coincide. 

Though  a  relatively  large  number  of  design  iterations  are 
required  for  convergence  in  Fig.  6,  the  improvements  in 
the  flutter  speed  X*  are  substantial.  This  objective  has  been 
increased  by  a  factor  of  3  over  the  baseline,  though  this  may 
not  be  a  logical  comparison,  as  the  baseline  is  composed  of 
a  fictitious  material  with  an  intermediate  density.  The  opti¬ 
mal  flutter  point,  474.08,  is  lower  than  the  unheated  flutter 
point  of  the  solid  panel  in  Fig.  3  (633.05),  though  this  isn’t 
a  good  comparison  either,  as  the  structure  in  Fig.  3  is  much 
heavier  (V*  =  1)  than  the  current  case  (V*  =  0.35).  The 
best  comparison  is  with  the  solid  panel  drawn  in  Fig.  4:  the 
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Fig.  5  Flutter-optimal  topology 
formation:  V*  =  0.35, 
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dynamic  eigenvalues  for  this  case  are  shown  in  Fig.  2,  with 
a  flutter  speed  of  40.52.  As  such,  the  optimal  topology  has  a 
flutter  speed  over  ten  times  higher  than  the  equivalent  solid 
panel  with  the  same  mass. 

For  further  comparison  with  Fig.  2,  the  eigenvalue  migra¬ 
tion  of  the  optimal  topology  is  given  in  Fig.  7.  The  flutter 
point  of  474.08  is  noted,  as  are  the  active  frequency  sepa¬ 
ration  constraints.  The  critical  separation  of  modes  2  and  3 
and  modes  5  and  6  occur  at  the  flutter  point,  the  critical  sep¬ 
aration  of  modes  4  and  5  occurs  just  before  the  flutter  point. 
For  the  solid  panel  of  Fig.  2,  the  higher  mode  eigenvalues 


experience  very  little  drift  as  A  increases,  the  imaginary  por¬ 
tion  of  modes  1  and  2  strongly  coalesce,  and  the  damping 
slope  of  the  flutter  mode  (dG/dX)  is  very  steep.  Contrast¬ 
ingly,  both  the  real  and  the  imaginary  portions  of  the  optimal 
eigenvalues  in  Fig.  7  are  strong  functions  of  A..  No  coales¬ 
cence  of  modes  1  and  2  occurs,  though  the  flutter  mode 
(i/rR*)  is  clearly  a  combination  of  these  two  modes  (shown 
in  Fig.  8),  where  larger  vibration  amplitudes  occur  near 
the  three-quarters  point  of  the  panel  length  (Dowell  1966). 
Furthermore,  the  damping  slope  at  flutter  is  very  shallow: 
if  the  panel  were  to  operate  at  a  flow  speed  just  beyond 


Fig.  6  Convergence  metrics 
of  the  topologies  in  Fig.  5: 
objective  function  A*  (top), 
frequency  separation 
constraints  (middle), 
and  damping  separation 
constraints  (bottom) 
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Fig.  7  Normalized  eigenvalue 
migration  for  the  optimal 
topology  of  Fig.  5 


*  flutter  point - eigenvalues  - - - active  Aa>  constraint  x  |q'3 


A.*,  unacceptable  deformations  would  accumulate  through¬ 
out  the  structure  relatively  slowly  (as  the  system  is  lightly 
unstable),  which  is  in  stark  contrast  to  most  aeronautical 
structures  (Bisplinghoff  et  al.  1955),  including  the  results  of 
Fig.  2. 

It  can  also  be  seen  in  Fig.  7  that  a  strong  coalescence  of 
modes  5  and  6  occurs  after  A.*:  mode  5  flutters  at  a  value 
of  521.12,  and  similarly,  mode  2  flutters  at  598.06.  These 
are  both  “conventional”  flutter  points  in  the  sense  that  the 
damping  slope  dG/dk  is  very  steep,  and  they  occur  via 
a  frequency  coalescence.  Without  the  presence  of  the  fre¬ 
quency  separation  constraints  in  (43),  these  higher  flutter 
modes  would  become  critical  during  the  optimization  pro¬ 
cess  (presumably  near  iteration  30,  which  is  when  the  first 
set  of  constraints  becomes  active  in  Fig.  6),  and  the  resulting 
design  space  discontinuity  would  hinder  further  improve¬ 
ments  in  A.*  (Odaka  and  Furuya  2005).  This  would  result  in 
a  substantially  non-optimal  topology,  as  A,*  is  increased  by 
over  100  between  iteration  30  and  the  final  design. 

The  presence  of  these  higher  mode  flutter  points  beyond 
A,*  =  474.08  are  of  no  consequence  in  a  deterministic 
sense,  as  the  flow  speeds  over  aerospace  structures  are  typ¬ 
ically  ramped  up  from  A.  =  0,  and  so  the  lowest  instability 
is  encountered  first  (Bisplinghoff  et  al.  1955).  In  a  non- 
deterministic  sense  however,  where  the  material  properties, 
topology  boundaries,  flow  details,  etc.,  may  be  inherently 
uncertain,  these  higher  flutter  modes  may  become  criti¬ 
cal  in  certain  areas  of  the  random  variable  space.  It  is 
demonstrated  by  Odaka  and  Furuya  (2005)  that  more  robust 
designs  are  obtained  with  higher  values  of  e<y,  via  smaller 
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Fig.  8  First  two  free- vibration  modes  and  the  flutter  mode  (irR*)  for 
the  optimal  topology  of  Fig.  5 


standard  deviations  in  the  flutter  speed.  Though  stochastic 
effects  are  not  considered  here,  the  topological  effect  of  e & 
is  demonstrated  in  Fig.  9,  where  the  topology/eigenvalues 
at  the  center  of  the  figure  are  reproduced  from  above. 
Increasing  the  required  separation  distance  from  0.95  to  a 
normalized  value  of  1.32  (i.e.,  making  the  constraint  harder 
to  satisfy)  decreases  the  flutter  speed  A.*  and  also  forces 
more  of  the  frequency  separation  constraints  to  become 
active,  as  expected.  The  next  flutter  point,  however,  is  well 
beyond  the  marked  value  of  A,*,  and  presumably  this  design 
(which  removes  material  from  the  center  mass  to  increase 
the  cross-bracing  at  the  edges  of  the  panel)  is  very  robust  in 
the  presence  of  uncertainties.  Opposite  trends  are  seen  for  a 
normalized  separation  distance  of  0.56,  which  has  the  high¬ 
est  flutter  speed,  the  fewest  number  of  active  constraints,  but 
a  modal  coalescence  (modes  3  and  4)  shortly  after  A.*. 

Fixing  the  normalized  frequency  separation  to  the  origi¬ 
nal  value  of  0.95,  the  effect  of  the  volume  fraction  V*  is  now 
explored.  The  upper  and  lower  areas  of  the  panel  which  are 
fixed  as  solid  account  for  a  volume  fraction  of  0.15,  so  the 
lowest  constraint  boundary  considered  in  Fig.  10  is  0.20.  For 
this  volume  the  topology  is  entirely  symmetric,  with  sparse 
truss-like  members  at  the  ends  and  center  of  the  panel. 
Increasing  the  volume  constraint  has  the  expected  result  of 
increasing  the  flutter  speed,  by  increasing  the  thickness  of 
the  cross-bracing  at  either  end,  as  well  as  the  width  of  the 
center  mass.  Topological  asymmetries  are  seen  for  every 
panel  above  V*  =  0.20  (where  the  design  for  V*  =  0.35  is 
repeated  from  above),  but  no  clear  pattern  emerges  concern¬ 
ing  which  side  of  the  panel  may  be  allocated  more  stiffness 
or  mass,  as  the  volume  is  increased.  The  final  topology  in 
Fig.  10  has  a  volume  of  0.587,  but  the  volume  constraint  for 
this  case,  which  was  set  to  0.60,  is  inactive.  Higher  values  of 
V*  therefore  need  not  be  considered,  as  the  optimizer  can¬ 
not  take  advantage  of  the  relaxed  constraint.  Similar  trends 
are  noted  in  thickness  design  variables  for  flutter-optimal 
panels:  areas  of  low  stiffness  increase  the  flutter  speed 
via  decreased  inertial  loads  (Weisshaar  1976;  Barboni  et  al. 
1999). 

This  point  is  reinforced  by  noting  that  all  topologies 
above  V*  =  0.40  have  a  higher  flutter  speed  than  the  fully 
solid  panel  of  Fig.  3  (A.*  =  633.05),  despite  weighing  much 
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Fig.  9  Flutter-optimal  *  flutter  point  - eigenvalues  - — ■ - active  dto  constraint 
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less.  The  flutter  speeds  of  solid  panels  of  volume  fractions 
between  0.05  and  1.00  (referencing  Fig.  4)  are  also  given  on 
the  right  of  Fig.  10;  each  of  which  is  non-Pareto-optimal  as 
compared  to  the  optimized  topologies.  The  nonlinear  rela¬ 
tionship  between  V*  and  X*  for  solid  panels  is  due  to  the  h? 
term  in  the  definition  of  the  bending  stiffness  D.  It  should 
be  noted  that  if  the  true  thickness  of  the  panel  V*  •  h  were 
used  in  the  definition  of  D ,  the  flutter  speed  X*  would  not 
be  a  function  of  V*,  instead  everywhere  equal  to  the  plot¬ 
ted  value  at  V*  =  1.00  (notwithstanding  finite  element 
discretization  errors).  It  can  also  be  seen  that  the  flutter 
speed  at  V*  =  1.00  has  a  slightly  different  flutter  speed 
than  that  quoted  for  Fig.  3.  This  is  due  to  the  fact  that  the 
data  in  Fig.  10  includes  structural  damping,  whereas  the 
verification  study  in  Fig.  3  does  not. 

The  topological  boundaries  in  most  of  the  panel  struc¬ 
tures  in  Fig.  10  are  well-defined,  though  “grey  stripes”  are 
present  in  three  of  the  cases:  V*  =  0.30,  0.40,  and  0.45. 
These  regions  of  intermediate-density  are  entirely  due  to  the 
filter-based  continuation  process.  When  the  radius  of  the 
filter  1-L  is  decreased  from  one  size  to  a  slightly  smaller 
size,  a  level  of  non-smooth  discreteness  is  introduced  to 
the  problem,  as  the  number  of  elements  which  fall  within 
the  neighborhood  of  any  central  element  drops.  This  effect 
causes  the  optimizer,  in  some  cases,  to  struggle  to  hold 
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the  highly  nonlinear  frequency  separation  constraints  A cop 
(many  of  which  are  active  during  the  continuation  process) 
directly  after  the  filter  radius  is  decreased,  and  some  grey 
material  is  added  in  response.  More  sophisticated  filters  and 
continuation  schemes  (which  gradually  force  the  topology 
to  entirely  black  and  white  designs  without  changing  the  fil¬ 
ter  size,  as  reviewed  by  Sigmund  2007)  may  out-perform 
the  methods  used  here,  particularly  with  regards  to  this 
problem. 

6.2  Thermal  buckling 

Next,  results  are  provided  for  buckling-optimal  topologies 
designed  under  (40),  where  aeroelastic  effects  are  not  con¬ 
sidered.  A  topological  convergence  for  a  volume  fraction 
V*  of  0.35  and  an  eigenvalue  separation  sp  of  0.5  is  seen 
in  Fig.  11.  For  the  initial  results  in  this  section,  the  heating 
is  entirely  uniform,  obtained  by  setting  Tp  =  T\j  —  1.  This 
effectively  turns  the  conduction  problem  (3)  “off’,  as  the 
temperature  of  each  finite  element  is  equal  to  unity,  regard¬ 
less  of  its  element  density.  All  terms  associated  with  dTx/dx 
in  the  sensitivity  analysis  derived  above  may  be  disregarded 
for  uniform  heating.  Comparing  Fig.  1 1  with  similar  con¬ 
vergence  results  in  Fig.  5,  it  can  be  seen  that  the  optimal 
buckling  problem  converges  much  quicker  than  the  optimal 


Fig.  10  Flutter-optimal 
topologies  for  increasing  volume 
fraction  (left),  and  the  Pareto 
front  (right):  e^/coo  =  0.95 
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Fig.  11  Buckling-optimal 
topology  formation  under 
uniform  heating:  V*  =  0.35, 
Sj  =  0.5 
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flutter  problem,  as  it  has  fewer  constraints.  Furthermore, 
these  eigenvalue  separation  constraints  (7*  —  7^_i  >  £7) 
are  only  active  for  higher  modes,  as  seen  in  Fig.  12.  It 
is  understood  that  the  normalized  buckling  eigenvalue  for 
mode  1  defines  the  normalized  objective  function  7?*,  and  at 
no  point  during  the  optimization  is  any  higher  mode  poised 
to  become  critical.  Modes  3  and  above  are  all  separated  by 
the  critical  distance  £7.  Numerical  experiments  (not  shown 
here)  indicate  that  increasing  77  by  an  order  of  magnitude 
has  no  discernible  effect  upon  the  topology.  This  lack  of 
dependence  on  eigenvalue  separation  is  in  stark  contrast  to 
the  flutter  results  (e.g.,  Fig.  9),  and  is  presumably  (from 
a  buckling  standpoint)  due  to  the  high  aspect  ratio  of  the 
geometry  considered  here  (Bendspe  and  Sigmund  2003). 

The  topology  in  Fig.  11  is  fully  converged  after  300 
iterations;  the  continuation  process  largely  removes  the 
grey  transition  areas  between  solid  and  void,  and  further 
increases  the  normalized  objective  function  7?*  as  well.  A 
relatively  high  normalized  value  of  7?*/^ 2  =  10.21  is 
obtained,  much  higher  than  the  value  of  4  quoted  above  for 
a  fully  solid  (and  much  heavier)  panel.  Objective  function 
improvements  over  the  baseline  design  are  moderate,  but  as 
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Fig.  12  Normalized  buckling  eigenvalues  corresponding  to  the 
topologies  in  Fig.  1 1 


above,  this  comparison  is  of  limited  usefulness  due  to  the 
use  of  intermediate-density  material.  The  optimal  topology 
is  entirely  symmetric  (about  the  mid-chord)  as  the  biasing 
effect  of  the  flow  velocity  is  not  present  here,  and  consists  of 
a  series  of  diagonally-oriented  truss-like  members.  Whereas 
the  flutter-optimal  panel  of  Fig.  5  utilizes  a  large  inertial 
mass  within  the  center  of  the  structure,  static-stress  metrics 
drive  the  design  process  in  Fig.  11.  According  to  the  critical 
buckling  mode  shown  in  Fig.  13,  flexural  stresses  will  be 
relatively  small  at  the  panel  center,  and  so  material  can  be 
removed  from  the  panel  center  and  reallocated  to  the  highly- 
stressed  panel  ends.  Largely  due  to  this  effect,  a  strong 
trade-off  between  buckling-  and  flutter-optimal  topologies 
is  to  be  expected,  and  will  be  demonstrated  in  the  next  sec¬ 
tion.  It  is  also  noted  that  the  third  buckling  mode  in  Fig.  13, 
though  largely  of  a  global  nature,  shows  some  local  buck¬ 
ling  patterns  along  the  thin  members  that  make  up  the  lower 
panel  surface. 

The  relationship  between  volume  fraction  and  the 
buckling-optimal  topology  is  seen  in  Fig.  14,  again  starting 
with  a  F*  of  0.20.  Like  the  flutter  data  of  Fig.  10,  beyond  a 
certain  threshold  (0.628  in  this  case)  the  volume  constraint 
boundary  becomes  inactive,  and  so  higher  values  are  of  no 
interest.  Indeed,  the  last  topology  in  Fig.  14,  which  was 
designed  under  V*  =  0.65,  is  non-Pareto-optimal  (in  the 
V*9  7?*  space)  with  respect  to  the  preceding  lighter  optimal 
topologies.  The  remaining  topologies  are  all  Pareto-optimal, 


mode  !:  Ryu  =  10.21 


mode  2;  =  I L73 


mode  3:  RyV  1 7.71 


Fig.  13  First  three  buckling  modes  of  the  optimal  topology  in  Fig.  1 1 
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Fig.  14  Buckling-optimal 
topologies  for  increasing 
volume  fraction  (left), 
and  the  Pareto  front  (right) 
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and  highly  superior  (for  a  given  volume  fraction)  to  the  crit¬ 
ical  buckling  temperature  of  the  solid  panels,  which  is  also 
plotted  in  Fig.  14.  As  before,  metrics  for  the  solid  panels 
are  obtained  using  concepts  of  Fig.  4.  Unlike  the  flutter 
data  of  Fig.  10,  the  buckling-optimal  topology  changes  in  a 
consistent  manner  with  increases  in  V*,  with  thickness  pri¬ 
marily  added  to  the  ends  and  the  center  of  the  panel  along 
with  a  continual  decrease  in  the  length  scale  of  the  diagonal 
truss-like  members. 

For  the  remainder  of  this  paper,  the  volume  fraction  is 
fixed  at  0.35.  Attention  is  now  turned  to  the  effect  of  a 
thermal  gradient;  leaving  the  upper  surface  temperature  Tjj 
at  unity,  the  buckling  topology  optimization  problem  of 
(40)  is  solved  for  various  values  of  the  lower  surface  tem¬ 
perature  Tl  between  zero  and  one.  Tl  values  above  one 
are  probably  unrealistic,  as  aerodynamic/radiated  heating 
will  typically  flow  from  the  environment  into  the  structure 
(Thornton  1992)  through  the  top  surface,  as  the  bottom  sur¬ 
face  is  not  exposed  to  the  flow.  It  should  also  be  noted 
here  that  radiation  heat  transfer  within  the  panel,  a  nonlin¬ 
ear  effect  which  may  be  important  at  elevated  temperatures, 
has  not  been  included  here.  Radiation  can  be  systemati¬ 
cally  included  in  topology  optimization  problems  with  a 
structure  surrounded  by  a  hot  environment  (a  cooling  fin, 
for  example)  via  a  nonlinear  convective  boundary  condi¬ 
tion  (Bruns  2007).  The  current  internal  radiation  case  is  far 
more  challenging  however,  with  the  role  of  radiation  view 
factors  in  the  variable-density  SIMP  approach  a  particularly 
interesting  complication. 


Topologies  for  Tl  values  of  0,  0.25,  0.5,  0.75,  and  1.0 
are  given  in  Fig.  15.  The  critical  value  of  R*  for  each 
design  is  given  as  well;  this  is  the  critical  buckling  mea¬ 
sured  at  the  design  temperature.  For  different  values  of  Tl, 
R*  for  a  given  design  will  change,  as  the  stresses  will  have 
been  redistributed  (4)  and  (7).  Typically,  decreasing  Tl  will 
increase  the  multiplicative  buckling  temperature  R*  due  to 
the  subsequent  decrease  in  stress  throughout  the  panel.  The 
upper  topology  in  Fig.  15  (7T  =  Tjj  =  1)  is  repeated 
from  above,  and  every  finite  element  has  the  same  temper¬ 
ature.  Lower  values  of  Tl  bring  about  minor  changes  in 
the  topology  at  the  center  and  the  ends  of  the  panel,  and 
the  number  of  thin  diagonal  truss-like  members  continually 
increases. 

Given  the  large  aspect  ratio  of  the  panel  and  the  spatially- 
uniform  upper  and  lower  surface  boundary  conditions,  the 
thermal  profile  for  most  cases  in  Fig.  15  is  entirely  expected: 
bulk  heat  transfer  from  top  to  bottom,  with  the  orientation 
of  the  temperature  gradient  vector  altered  to  travel  along 
the  local  axis  of  individual  truss-like  members.  In  light  of 
this,  an  interesting  question  is  the  impact  that  simply  fix¬ 
ing  the  temperature  distribution  (a  linearly-varying  spatial 
distribution  between  Tl  and  Tu,  for  example)  would  have 
on  the  topology  optimization.  Taking  the  case  of  Tl  =  0.5 
as  a  representative  example,  setting  X  j  equal  to  zero  in 
(30)  (which  essentially  removes  the  relationship  between 
the  elemental  density  and  the  temperature  distribution)  will 
change  dTk/dx  by  7.2  %  for  the  initial  baseline  design,  and 
12.3  %  for  the  final  optimal  design  in  Fig.  15:  a  moderate 


Fig.  15  Buckling-optimal 
topologies  and  temperature 
fields  for  increasing  values  of  Tl 
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impact.  This  result  is  problem-dependent  of  course,  depend¬ 
ing  somewhat  upon  the  geometry  of  the  design  domain,  and 
strongly  upon  the  type  of  boundary  condition  (convection 
heat- flux,  e.g.). 

For  a  Tl  of  zero  however  (and  to  a  lesser  extent,  0.25), 
the  optimizer  has  developed  a  topology  in  which  the  upper 
heated  surface  is  entirely  disconnected  from  the  remainder 
of  the  panel,  which  thus  remains  unheated  and  unstressed. 
The  critical  eigenvalue  for  this  case  is  purportedly  very 
high  (R* /tv2  =  19.22),  though  clearly  the  buckling  eigen¬ 
value  for  the  upper  heated  layer,  considered  on  its  own, 
should  be  much  lower  (i.e.,  the  buckling  mode  should  be 
entirely  localized  to  the  upper  layer,  but  is  instead  of  a  global 
nature).  A  closer  inspection  reveals  that  the  upper  surface 
and  the  remainder  of  the  truss  members  are  separated  by  a 
single  row  of  low  density  ( xe )  finite  elements.  Their  den¬ 
sity  is  low  enough  such  that  very  little  thermal  energy  may 
pass  through,  but  high  enough  to  structurally  reinforce  the 
surface  at  several  locations.  Clearly,  this  topology  is  of  lit¬ 
tle  practical  value;  the  optimizer  has  taken  advantage  of 
the  problem  description  in  general  and  the  finite  element 
discretization  in  particular.  Additional  constraints  may  be 
formulated  into  (41)  to  prevent  this  behavior,  but  as  will  be 
seen  in  the  next  section,  this  structure  is  highly  non-optimal 
in  an  aerothermoelastic  sense. 

6.3  Flutter  of  a  heated  panel 

Flutter-optimal  results  are  now  provided  for  heated  panels, 
which  involves  solving  the  optimization  problem  of  (43) 
for  various  value  of  the  multiplicative  temperature  T .  It 
is  repeated  here  that  T  is  in  no  way  related  to  the  buck¬ 
ling  temperature  T*  of  the  panel’s  topology  (which  is  not 
included  in  the  design  optimization),  but  instead  is  a  quan¬ 
tity  which  is  held  fixed  during  the  optimization.  A  series  of 
topologies  are  given  in  Fig.  16  for  various  levels  of  uniform 
heating,  along  with  the  value  of  R  used  to  design  the  struc¬ 
ture,  the  resulting  objective  function  (flutter  speed,  A,*),  and 
the  resulting  critical  buckling  load  7?*.  This  latter  metric  is 
not  utilized  by  the  optimizer,  but  is  provided  for  complete¬ 
ness.  The  first  topology  in  Fig.  16  is  repeated  from  Fig.  5, 
which  was  developed  under  zero  thermal  load.  The  buckling 
temperature  of  this  case  is  R* /tv2  =  6.03,  which  is  much 
lower  than  the  buckling-optimum  result’s  in  Fig.  11.  This 


would  confirm  the  topological  trade-off  between  buckling- 
and  flutter-optimal  structures  discussed  above. 

Increasing  the  value  of  the  applied  temperature  R  in 
Fig.  16  has  a  substantial  effect  upon  the  flutter-optimal 
topology.  With  increasing  R ,  the  inertial  mass  at  the  cen¬ 
ter  of  the  panel  becomes  less  defined  (and  is  completely 
removed  from  the  hottest  structures),  and  the  cross-bracing 
is  increased  at  the  ends  of  the  panel.  The  maximum  attain¬ 
able  flutter  speed  generally  decreases,  as  higher  thermal 
loads  weaken  the  structure,  but  the  critical  buckling  load 
R* /tv2  is  relatively  insensitive,  between  6.0  and  6.4  for 
each  case.  This  is  surprising  in  the  sense  that  the  topologies 
optimized  for  flutter  under  higher  thermal  loads  clearly  are 
being  driven  by  strategies  observed  for  buckling-optimum 
structures  in  the  previous  section.  The  third  structure  of 
Fig.  16  is  an  outlier  in  the  sense  that  its  flutter  point  is 
larger  than  the  adjacent  structures  optimized  for  higher  or 
lower  heating  levels.  This  may  be  indicative  of  a  highly  - 
complex  design  space,  or  the  fact  that  the  continuation 
process  in  this  case  has  not  provided  a  truly  solid-void 
design.  The  optimizer  may  be  taking  advantage  of  this  fic¬ 
tional  intermediate-density  material  (seen  in  the  center  of 
the  panel)  to  obtain  a  higher  flutter  speed  than  would  be 
obtained  if  every  finite  element  was  explicitly  forced  to  its 
density  limits. 

Off-design  behavior  for  the  five  topologies  of  Fig.  16,  the 
buckling-optimum  design  in  Fig.  11,  and  a  solid  panel,  are 
given  in  Fig.  17.  Each  panel  in  this  figure  has  the  same  vol¬ 
ume  (V*  =  0.35),  and  each  is  largely  composed  of  solid 
and  void  elements.  The  plot  is  provided  in  the  same  format 
as  the  verification  study  of  Fig.  3:  the  flutter  boundary  A,* 
is  given  has  a  function  of  the  applied  thermal  load  R ,  and 
the  boundary  across  which,  for  a  given  value  of  R ,  the  panel 
is  blown  flat,  is  also  shown.  This  latter  boundary  emanates 
from  the  buckling  load  7?*,  and  eventually  merges  with  the 
flutter  boundary  at  some  larger  value  of  R.  The  design  point 
in  (7?,  A.*)  space  is  also  indicated  in  the  figure.  Several  inter¬ 
esting  observations  can  be  made  from  Fig.  17.  First,  many 
topologies  are  seemingly  out-performed  by  another  topol¬ 
ogy’s  off-design  behavior.  For  example,  the  unheated  panel 
has  a  flutter  speed  of  474.08  at  its  design  point  (R/tv2  =0), 
but  the  panel  optimized  at  R/tv2  =  6.11  has  a  higher  flutter 
speed  (541 .93)  at  this  same  point.  Only  the  unheated  panel  is 
feasible  at  this  point  however,  as  the  other  topologies  violate 


Fig.  16  Flutter-optimal 
topologies  under  various  levels 
of  uniform  heating:  V*  =  0.35, 
£cd/m0  =  0.95 
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Fig.  17  Off-design  flutter  and 
buckling  boundaries  for  panels 
with  a  volume  fraction  of  0.35, 
under  uniform  heating 
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the  frequency  separation  constraints  when  their  off-design 
behavior  is  computed,  and  as  such,  cannot  be  considered  as 
robust. 

An  indicator  of  this  robustness  manifests  itself  through 
changes  in  the  critical  mode  along  the  flutter  boundary.  At 
the  design  point,  the  critical  flutter  mode  is  always  the  first 
mode,  as  this  is  stipulated  by  the  frequency  and  damping 
separation  constraints  added  to  the  optimization  problem 
(43).  For  an  off-design  point,  however,  the  critical  flutter 
mode  could  be  due  to  an  interaction  of  higher  modes.  Again 
taking  the  panel  designed  under  unheated  conditions  as  an 
example:  the  flutter  speed  starts  at  the  designed  value  of 
474.08,  and  decreases  in  a  monotonic  fashion  with  increas¬ 
ing  values  of  R ,  which  is  very  similar  to  the  behavior  seen 
for  a  solid  panel  (Fig.  3).  At  R/n2  =  3.53,  a  coales¬ 
cence  of  modes  2  and  3  becomes  more  critical  than  the 
original  flutter  mode,  and  A*  drops  off  at  a  much  faster 
rate.  At  R/n2  =  4.91,  the  flutter  mode  changes  again 
(a  coalescence  of  modes  5  and  6),  after  which  the  flutter 
speed  actually  increases  with  increased  heating.  This  would 
indicate  a  highly  complex  relationship  between  the  flutter¬ 
ing  fifth  vibration  mode  and  the  stresses  accumulated  by 
uniform  heating. 

At  R /n 2  =  6.07,  the  flutter  mode  reverts  back  to  the  orig¬ 
inal  first  mode,  and  X*  decreases  until  the  flutter  boundary 


merges  with  the  buckling  boundary.  The  end  result  is  a  dan¬ 
gerous  “flutter  dip”,  whose  presence  would  go  undetected 
if  higher  mode  effects  were  not  included  in  the  analy¬ 
sis.  A  discontinuous  flutter  boundary  (namely,  a  dip)  via 
mode- switching  is  a  commonly-seen  behavior  in  aeroelas- 
tic  systems  undergoing  parametric  sweeps:  see  for  example 
Forster  and  Yang  (1998)  (via  changes  in  skin  thickness)  and 
Beran  et  al.  (2004)  (via  changes  in  Mach  number),  among 
many  others. 

In  many  of  the  cases  in  Fig.  17,  these  areas  of  low  flut¬ 
ter  occur  only  for  highly  off-design  areas  of  the  plot.  An 
alternative  aerothermoelastic  optimization  procedure  to  the 
one  taken  here  is  to  maximize  the  unheated  flutter  speed 
X*  under  a  series  of  constraints  upon  the  critical  buckling 
temperature  7?*:  in  this  way,  a  Pareto  front  may  be  formed. 
The  aeroelastic  behavior  of  these  designs  under  (off-design) 
heated  conditions,  however,  is  likely  to  be  very  poor,  due  to 
the  advent  of  higher  flutter  modes  and  their  associated  dips. 
The  procedure  used  here,  to  maximize  flutter  directly  under 
a  heated  environment,  is  superior.  The  panel  optimized  at 
R/n2  =  6.11,  for  example  (bottom  topology  of  Fig.  16), 
has  a  relatively  high  flutter  speed  through  the  entire  heat¬ 
ing  range  of  Fig.  17,  and  a  very  mild  dip  between  1.21  and 
4.37.  As  expected,  the  flutter  performance  of  the  buckling- 
optimal  design  is  poor  despite  having  the  largest  7?*  in  the 


Fig.  18  Normalized  eigenvalue 
migration  for  the  topology 
optimized  at  R/n2  =6.11: 
increased  X  for  a  fixed  heating 
of  R/n2  =  7.5 
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Fig.  19  Flutter-optimal 
topologies  under  various  levels 
of  non-uniform  heating: 

V*  =  0.35,  e^jcDo  =  0.95 
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plot,  and  both  the  heated  flutter  and  the  buckling  character¬ 
istics  of  the  solid  panel  (with  a  volume  fraction  of  0.35)  are 
very  low. 

The  qualitative  shape  of  each  design’s  lower  buckling 
boundary  is  the  same,  which  would  indicate  that  higher 
buckling  modes  do  not  become  active  in  the  same  way  that 
is  observed  along  the  flutter  boundary.  To  further  explore  the 
“blown  flat”  concept  (which  occurs  across  this  boundary), 
eigenvalue  migration  for  the  flutter-optimal  panel  topology 
designed  at  R/tt2  =  6.11  are  given  in  Fig.  18.  This  is  done 
across  a  range  of  flow  speeds  A,  for  a  fixed  value  of  the 
heating  parameter  R/n2  =  7.50,  which  is  greater  than  the 
buckling  temperature  of  this  design  (R*/jr2  =  6.36,  given 
in  Fig.  16).  For  low  values  of  k  the  panel  is  buckled:  one 
of  the  eigenvalues  has  a  positive  real  part,  but  its  imaginary 
portion  (frequency)  is  0,  and  so  this  is  a  static  instability. 
The  panel  is  blown  flat  at  k  =  37.50,  when  the  imaginary 
portion  becomes  positive  and  the  real  part  negative  (stabi¬ 
lizing)  This  occurs  due  to  the  collision  of  two  real-valued 
eigenvalues,  which  then  become  complex  conjugates.  This 
same  pair  loses  stability  at  k*  =  175.68,  which  defines  the 
flutter  speed. 

The  A. -marching  concept  outlined  above  to  locate  a  good 
initial  guess  for  the  flutter  point  will  fail  in  the  case  of 
Fig.  18.  To  repeat,  k  is  increased  from  a  low  value  (typ¬ 
ically  0)  until  G  =  max(Re(Pk ))  becomes  positive.  But 
it  can  be  seen  in  the  figure  that  G  is  positive  at  k  =  0 
(due  to  the  static  buckling  instability),  though  this  clearly 
is  not  the  desired  flutter  point.  The  problem  is  solved  by 


only  considering  eigenvalues  whose  imaginary  portions  are 
greater  than  zero  (i.e.,  dynamic)  in  the  identification  of  G. 
This  issue  is  not  encountered  for  any  of  the  cases  given  in 
Fig.  16,  as  none  are  operating  under  a  thermal  load  greater 
than  their  buckling  load,  and  so  all  of  their  eigenvalues  will 
have  a  nonzero  imaginary  part.  The  exercise  of  Figs.  16  and 
17  is  repeated  for  the  case  of  non-uniform  heating.  Two  of 
the  topologies  are  subjected  to  a  temperature  in  the  post- 
buckled  regime,  and  so  the  modification  to  the  A. -marching 
process  is  required.  Flutter-optimal  topologies  under  non- 
uniform  heating  (specifically,  Tl  =  0)  are  given  in  Fig.  19 
and  off-design  behavior  is  given  in  Fig.  20. 

The  first  topology  in  Fig.  19  is  repeated  from  above, 
designed  without  a  thermal  load.  The  critical  buckling  tem¬ 
perature  is  higher  than  that  cited  in  Fig.  16,  as  less  of  the 
panel  is  being  heated,  though  of  course  the  flutter  speed  is 
the  same.  As  above,  increasing  the  thermal  load  on  the  struc¬ 
ture  leads  the  optimizer  to  remove  the  central  mass  in  favor 
of  additional  cross-bracing  at  the  panel  ends.  For  the  high¬ 
est  thermal  load,  the  structure  is  almost  entirely  symmetric, 
composed  of  a  series  of  thin  diagonally-oriented  members. 
Despite  the  fact  that  flutter  speed  is  the  actual  objec¬ 
tive  function,  the  buckling- sensitive  aspects  of  the  problem 
drive  the  design  process:  in  (18),  the  energy  associated 
with  the  thermal  term  T  ■  Ka  dominates  the  aerodynamic 
contribution  from  k-  Ka  during  the  A.*  computation.  A  com¬ 
promise  is  forced  between  flutter-optimal  topologies  and 
bucking-optimal  topologies,  although  the  latter  metric  is  not 
explicitly  included  in  the  optimization  statement. 


Fig.  20  Off-design  flutter  and 
buckling  boundaries  for  panels 
with  a  volume  fraction  of  0.35, 
under  non-uniform  heating 


*  design  point 

optimized  ai  R/tt2  =  0.00 

- optimized  ai  RJn2  =  9.17 

- optimized  at  R/rc2  =  1 2.22 

optimized  at  R/rc2  =  15.28 
- -  optimized  for  buckling 
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The  lower  surface  temperature  Tl  is  set  to  zero  for 
these  results,  a  potentially  poor  choice  given  that  the  buck¬ 
ling  optimizer  gave  questionable  results  for  this  value  in 
Fig.  15.  It  is  seen  in  Fig.  20  that  the  off-design  aerother¬ 
moelastic  performance  of  the  buckling-optimal  design  is 
very  poor  however,  with  consistently  low  flutter  speeds. 
Furthermore,  none  of  the  topologies  in  Fig.  19  mimic  the 
structures  seen  above,  where  the  heated  upper  panel  was 
largely  disconnected  from  the  remainder  of  the  topology. 
The  off-design  flutter  boundaries  in  Fig.  20  are  similar  to 
those  found  for  uniformly-heated  panels,  with  strong  flut¬ 
ter  dips  for  heating  values  far  less  than  the  design  point. 
It  is  further  noted  that  data  for  a  solid  panel  has  not  been 
included  in  Fig.  20:  because  the  lower  portion  of  this  panel 
is  fixed  as  entirely  void,  the  imposition  of  a  thermal  gra¬ 
dient  from  the  top  to  the  bottom  surface  is  not  a  feasible 
situation. 


7  Conclusions 

This  work  has  optimized  the  internal  topology  of  a  metallic 
panel  subject  to  a  high  supersonic  flow  over  its  upper  sur¬ 
face,  and  prescribed  temperature  boundary  conditions  along 
the  upper  and  lower  surfaces.  The  cross-section  of  the  panel 
is  discretized  into  finite  elements,  each  of  which  is  assigned 
a  density-based  design  variable  which  indicates  a  solid  ele¬ 
ment  or  a  void.  The  quasi-steady  temperature  distribution 
throughout  the  panel  is  computed,  followed  by  the  ther¬ 
mal  stresses  and  the  critical  buckling  temperatures.  Using 
well-known  piston  theory  aerodynamics,  the  flutter  speed 
of  the  structure  can  be  found  for  a  given  thermal  load, 
which  occurs  at  a  Hopf-bifurcation  point.  Both  eigenvalue- 
related  metrics  (buckling  temperature  and  flutter  point) 
are  optimized  using  the  SIMP-based  topology  optimiza¬ 
tion  methodology,  though  each  problem  must  be  outfitted 
with  highly  nonlinear  eigenvalue- separation  constraints  to 
prevent  switching  of  the  critical  mode  (and  the  associ¬ 
ated  design  space  discontinuities)  during  the  optimization 
process.  Following  a  model  verification  study,  three  sets 
of  results  are  presented:  unheated  flutter-optimal  topolo¬ 
gies,  thermal  buckling-optimal  topologies  (no  aerodynamic 
loading),  and  flutter  optimal  topologies  under  a  thermal 
load  (aerothermoelastic).  The  following  conclusions  can  be 
drawn: 

1.  Unheated  flutter-optimal  topologies  can  be  obtained 
by  allocating  mass  to  the  center  of  the  panel  (inertial 
effects)  and  cross-bracing  at  the  ends  of  the  panel  (elas¬ 
tic  effects).  Due  to  the  biasing  effect  of  the  flow  vector, 
asymmetries  develop  in  the  structure  later  in  the  opti¬ 
mization  process,  though  the  initial  iterates  are  entirely 
symmetric. 


2.  Whereas  the  flutter  of  a  solid  panel  occurs  very  rapidly 
(steep  damping  slope)  after  the  coalescence  of  the  imag¬ 
inary  portion  of  two  eigenvalues,  the  flutter  of  an  opti¬ 
mized  topology  does  not  involve  a  strong  coalescence, 
and  the  damping  slope  through  the  Hopf-bifurcation  is 
very  moderate. 

3.  Frequency  separation  constraints  are  typically  active  at 
values  of  k  both  equal  to  and  less  than  the  flutter  point. 
Increasing  the  required  separation  distance  drops  the 
optimal  flutter  value  and  increases  the  number  of  active 
constraints,  but  should  provide  more  robust  designs  in 
the  presence  of  uncertainties,  as  higher  mode  flutter 
points  are  forced  to  reside  at  much  higher  A, -values 
than  A,*. 

4.  The  thermal  buckling  optimization  problem  converges 
much  faster  than  the  flutter  problem,  as  fewer  eigen¬ 
value  separation  constraints  are  needed,  and  are  typ¬ 
ically  only  active  among  the  higher  buckling  modes. 
These  designs  involve  removing  all  mass  from  the  cen¬ 
ter  of  the  panel  (where  the  stresses  will  be  small) 
in  favor  of  truss-like  members.  These  topologies  are 
entirely  symmetric. 

5.  Accuracy  issues  are  noted  in  the  buckling-optimal 
topology  when  designed  under  a  large  thermal  gradient, 
as  the  optimizer  obtains  a  structure  where  the  heated 
upper  panel  is  entirely  disconnected  from  the  lower  por¬ 
tion  of  the  panel,  which  thus  remains  unheated.  This 
structure  is  found  to  be  highly  non-optimal  in  a  flutter 
sense. 

6.  Aerothermoelastic  conditions  are  obtained  by  optimiz¬ 
ing  the  flutter  speed  under  a  specified  thermal  load. 
Increasing  this  load  forces  a  compromise  between  the 
flutter-optimal  topologies  and  buckling-optimal  topolo¬ 
gies,  although  buckling  is  not  explicitly  included  as  an 
objective  function  or  a  constraint. 

7.  Off-design  flutter  behavior  for  these  structures  (over  a 
range  of  heating  levels)  is  characterized  by  changes  in 
the  critical  flutter  mode,  which  can  cause  drastic  drops 
in  the  flutter  speed.  Robust  panels  will  clearly  need  to 
simultaneously  consider  several  heating  conditions  dur¬ 
ing  the  topology  optimization  process,  a  common  tactic 
in  aerospace  design. 
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